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1.0  Introduction 


Considerable  progress  has  been  made  in  the  understanding  of  magnetospheric 
processes  in  the  last  25  years.  During  this  time  much  of  the  effort  has  been  focused 
on  understanding  processes  operating  in  the  tail  of  the  magnetosphere  and  near  the 
magnetospheric  boundaries.  The  inner  magnetosphere  has  not  been  extensively 
investigated  in  the  last  10  to  20  years.  The  Combined  Release  and  Radiation 
Effects  Satellite  (CRRES)  Program  is  designed  to  make  a  substantial  contribution  to 
understanding  this  region  of  space. 

In  this  effort  McDonnell  Douglas  Space  Systems  Company  (MDSSC)  will  introduce 
novel  new  modelling  approaches  and  will  satisfy  one  of  main  goals  of  the  CRRES 
mission,  the  development  of  new  static  and  dynamic  radiation  belt  models.  Such 
models  are  an  important  tool  for  engineers  for  the  design  of  systems  that  can 
survive  in  the  space  environment.  The  present  radiation  belt  models  that  are  now 
used  by  both  the  scientific  and  engineering  communities  are  the  Vette  radiation 
models  developed  by  the  National  Space  Sciences  Data  Center  (NSSDC)  in  the  late 
1 960’s  and  early  70's.  The  data  sets  which  were  used  to  develop  the  Vette  models 
were  acquired  by  instruments  that  are  quite  primitive  compared  to  todavs  state-of- 
the-art  instruments.  Nevertheless  these  older  models  have  served  the  community 
well. 

The  present  NSSDC  developed  models  are  organized  in  R,  L  space,  a  coordinate 
system  developed  by  C.  E.  Mcllwain  in  1961.  This  coordinate  system  has  been 
virtually  unchanged  since  that  time.  Improvements  have  come  only  in  the  form  of 
improved  computational  techniques.  Although  the  B,  L  system  has  proven  useful  in 
the  inner  zone,  its  use  in  the  outer  zone  has  not  been  as  successful. 
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This  CRRES  analysis  effort  will  include  the  development  of  new  tools  for  the 
organization  of  the  data.  This  effort  will  develop  novel  new  techniques  for 
organizing  charged  particle  data  in  the  inner  and  outer  zone.  In  the  inner  zone  we 
will  develop  a  model  that  not  only  takes  into  account  the  effect  of  the  magnetic  field 
in  organizing  the  charged  particles  but  also  the  effect  of  the  solar  cycle  dependent 
atmosphere  in  shaping  the  low  altitude  region  of  the  inner  radiation  belt.  In  the 
outer  zone  we  will  provide  a  coordinate  system  that  can  correctly  represents 
adiabatic  changes  in  the  radiation  belt  and  fully  takes  into  account  drift  shell  splitting 
and  yet  represents  the  entire  outer  zone  in  terms  of  only  two  parameters,  the  first 
and  second  invariant  for  each  observation  and  pitch  angle.  Once  the  new  tools 
developed  under  this  effort  are  implemented  and  a  best  fit  is  made  to  the  CRRES 
data,  a  high  quality  radiation  belt  model  will  be  created  that  will  be  valid  for  all  epoch 
within  a  solar  cycle  and  for  all  magnetic  conditions  of  the  magnetosphere.  In  the 
inner  zone,  the  new  model  will  permit  calculation  of  the  fluxes  during  any  part  of  the 
solar  cycle.  In  the  outer  zone  the  model  will  help  separate  adiabatic  changes  from 
non  adiabatic  variations,  allow  the  calculation  of  particle  fluxes  for  all  states  of 
magnetospheric  compression,  and  will  help  theorists  explain  many  of  the  observed 
changes  within  the  magnetosphere. 

The  analysis  effort  consists  of  two  distinct  phases,  the  first  of  these  is  the 
development  of  the  required  computer  code  for  defining  the  new  coordinate  system 
and  the  second  is  the  fitting  of  the  data  utilizing  the  new  coordinate  system.  By  its 
very  nature  the  development  of  the  computer  code  cannot  be  completely  decoupled 
from  the  analysis  of  the  data.  The  quality  of  the  data  and  the  various  features  found 
within  the  data  stream  will  dictate  the  ultimate  development  of  the  model  coordinate 
system  and  the  final  fit  of  the  data. 
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2.0  A  Pitch  Angle  Dependent  Invariant  Routine 


During  the  first  year  of  the  contract  almost  all  of  the  effort  has  been  spent 
developing  and  polishing  the  various  software  tools  that  will  be  required  to  produce 
a  new  high  precision  CRRES  model  of  the  inner  and  outer  zones.  Some  of  this 
software  must  yet  be  verified  using  CRRES  data.  This  technical  report  will  describe 
the  software  that  is  considered  fully  operational  and  which  is  available  for  use  by 
any  of  the  CRRES  investigators. 

This  annual  report  will  concentrate  on  describing  the  software  that  could  be  verified 
without  using  the  CRRES  data.  The  major  achievement  during  the  first  year  was 
the  development  of  a  new  B,L  code  that  calculates  the  invariant  at  the  satellite  for 
many  different  pitch  angle  particles.  Calculating  B,L  has  always  been  a 
computationally  expensive  procedure  and  thus  considerable  effort  was  expended  to 
produce  a  code  that  is  efficient  and  cost  effective.  The  calculation  of  the  second 
invariant  requires  the  calculation  of  a  line  integral  that  makes  many  calls  to  the 
magnetic  field  subroutine. 

2. 1  A  Fast  Version  of  IGRF 

Present  versions  of  B,L  use  only  the  internal  model  of  the  magnetic  field.  The 
CRRES  code  must  use  both  internal  and  external  models  of  the  magnetic  field, 
since  it  must  take  into  account  drift  shell  splitting  in  the  outer  magnetosphere.  Thus 
one  of  the  most  important  routines  for  saving  computer  time  is  the  development  of 
an  internal  and  external  magnetic  field  routine  that  optimizes  computer  speed.  The 
Olson-Pfitzer  1 977  tilt  dependent  model  is  such  a  routine.  It,  however,  uses  the 
Barraclough  internal  field  routine  and  thus  is  not  appropriate  for  the  CRRES  effort. 
The  IGRF  routines  using  the  modern  field  coefficients  were  obtained  from  the 
National  Space  Sciences  Data  Center  (NSSDC).  These  routines  are,  however, 
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considerably  slower  than  the  routine  contained  in  the  original  version  of  the  1 977  tilt 
dependent  model.  The  main  field  routine  contained  in  the  1977  model  is  derived 
from  Joe  Cain's  SPHRC  routine.  This  routine  gains  additional  speed  at  the  expense 
of  some  memory.  Instead  of  using  indexed  loops  it  explicitly  writes  out  the  spherical 
harmonic  expansion  terms  and  thus  all  of  the  overhead  required  to  keep  track  of  the 
various  indices  is  abolished.  It  is  this  authors  opinion  that  this  version  of 
representing  the  main  fieid  ts  inherently  50%  faster  than  any  other  representation. 
The  version  used  in  the  1977  tilt  dependent  model  has  the  Gauss  normalized 
Barraclough  coefficients  built  into  the  model.  It,  furthermore,  has  a  term  dropping 
algorithm,  that  drops  the  higher  order  terms  as  distances  increases.  This  results  in 
a  considerable  savings  in  computer  time. 

The  IGRF  internal  field  model  developed  for  the  CRRES  analysis  begins  with  Joe 
Cain's  spherical  harmonic  expansion.  The  new  IGRF  routine  has  been  given  the 
name  SPIGRF  (SPeed  IGRF).  Since  the  IGRF  coefficients  are  for  a  tenth  order 
expansion,  terms  up  the  N=1 1  are  contained  in  SPIGRF  (the  Barraclough  model 
only  had  10  terms).  The  first  time  SPIGRF  is  called,  it  calls  a  routine  called 
FLDCOF.  FLDCOF  reads  the  appropriate  IGRF  coefficient  sets,  interpolates  to  the 
epoch  of  interest  and  converts  the  Schmitt  normalized  IGRF  coefficients  to  Gauss 
normalized  coefficients.  If  during  a  subsequent  call,  the  date  has  changed  by  more 
than  0.1  year,  FLDCOF  is  once  again  called  to  update  the  coefficients  to  the  new 
epoch.  It  is  not  computationally  efficient  to  update  the  coefficients  for  minor 
changes  in  time.  In  fact  it  might  be  appropriate  for  the  CRRES  mission  to  select  a 
specific  date  and  use  coefficients  for  the  interna!  magnetic  field  that  do  not  change 
with  time.  This  is  a  programmatic  decision  and  will  have  little  or  no  impact  on  the 
science  unless  CRRES  continues  to  function  for  more  than  4  or  5  years. 
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The  term  dropping  algorithm  that  is  a  part  of  SPIGRF  is  a  smooth  algorithm.  The 
algorithm  uses  predetermined  altitudes  where  specific  terms  are  to  be  discontinued. 
Table  1  lists  the  altitudes  and  maximum  number  of  terms  that  are  used  by  SPIGRF. 

Table  1 


Altitude,  Re 

Number  of  terms  jj 

12.0 

2 

8.0 

3 

6.0 

4 

5.0 

5 

4.0 

6 

3.2 

7 

2.5 

8  I 

2.0 

9  1 

1.6 

10 

1.4 

11 

Thus  for  altitudes  greater  than  1 2  Re  only  the  dipole  term  is  used.  For  altitudes  less 
than  1 .4  Re  all  1 1  terms  are  used.  For  altitudes  between  1 .4  and  1 .6,  the 
contribution  of  the  N=1 1  term  is  linearly  reduced  from  its  full  value  at  1 .4  to  zero  at 
1 .6.  Similarly  for  all  other  intervals.  The  altitude  values  at  which  the  term  dropping 
takes  place  was  experimentally  determined  by  comparing  the  truncated  model  with 
the  untruncated  version.  The  truncated  model  differs  from  the  untruncated  model 
by  no  more  than  0.1  nanotesla. 
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This  new  IGRF  code  is  inherently  1 .5  times  as  fast  as  the  NSSDC  code,  FELDG, 
when  used  with  all  the  terms.  The  term  drop  off  algorithm  significantly  improves  the 
speed  without  sacrificing  accuracy.  Figure  1  presents  the  speed  advantage  of 
SPIGRF  when  compared  to  the  original  NSSDC  code.  At  3.0  Re  it  is  3  times  as  fast 
as  FELDG,  7  times  as  fast  at  5  Re,  10  times  as  fast  at  7  Re,  and  35  times  as  fast  at 
12  Re.  Since  the  term  drop  off  algorithm  removes  the  effect  of  a  term  smoothly, 
there  are  no  discontinuities  in  the  field.  For  a  typical  CRRES  orbit,  this  version  of 
IGRF  should  have  an  average  speed  advantage  of  7  or  8. 


Figure  1 

Improvements  in  IGRF  Calculations 

by  using  SPIGRF  with  term  dropping  algorithm 


0  2.5  5.0  7.5  10.0  4  ?  5 

Distance  in  Earth  Radi; 


A  listing  and  a  brief  description  of  the  calling  sequence  for  SPIGRF  is  given  in 

Appendix  A. 
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2.2  Internal  plus  External  Field 


The  new  routine,  SPIGRF,  was  combined  with  the  remaining  routines  of  the  1 977  tilt 
dependent  magnetic  field  model  to  produce  a  high  speed  quiet  time  magnetic  field 
model  that  uses  the  IGRF  coefficients.  The  name  of  the  entire  file  that  contains  the 
tilt  dependent  model,  SPIGRF,  and  a  test  routine  was  given  the  name  BMNIGRF. 

The  execution  speed  of  the  new  fast  IGRF  code,  SPIGRF,  plus  the  1977  tilt 
dependent  external  model  is  faster  everywhere  than  the  old  internal  IGRF  code, 
FELDG,  without  the  external  model.  At  1 .4  Re  or  less  the  speed  advantage  is  1 .2, 
at  2  Re  it  is  1 .05,  at  6  Re  it  is  1 .6  and  at  10  Re  it  is  2.0.  Thus  calculating  B,L  using 
the  external  and  internal  field  routines  is  faster  than  calculating  a  B,L  based  on  the 
internal  field  alone  using  the  older  IGRF  routines. 

Note:  The  external  field  model  presently  used  for  the  initial  CRRES  studies  is  the 
1977  quiet  time  tilt  dependent  model.  The  dynamic  model,  which  will  be  used  to 
study  the  outer  zone  as  a  function  of  magnetospheric  compression,  will  replace  the 
quiet  time  model  at  the  appropriate  time. 

Appendix  B  gives  a  listing  of  the  BMNIGRF  test  routine,  the  1977  tilt  dependent 
routine,  BXYZMU,  and  the  various  routines  required  to  combine  the  external  and 
internal  magnetic  fields.  The  internal  field  routines  are  described  in  Appendix  A. 

2.3  A  Pitch  Angle  Dependent  B,L  Routine 

In  order  to  adequately  represent  directional  data,  it  is  necessary  to  define  the  first 
and  second  invariant  for  each  directional  measurement.  At  the  present  time  the 
current  B,L  routines  calculate  the  invariant  for  a  particle  that  mirrors  at  the  location 
of  the  satellite.  All  other  particles  will  of  course  have  different  first  and  second 


invariant.  The  first  invariant  is  simply  the  mirror  point  magnetic  field  of  the  particle 
and  thus  Bmjr  is  given  by 

R  -  ^  local 

Dmir  ~  .  1.  . 

sin  (alocal)  (1) 

where  B|0ca|  is  the  magnetic  field  a:  the  location  of  the  measurement  and  otjoca|  is 

the  local  pitch  angle  of  the  measurement.  Thus  the  first  invariant  can  be  easily 
calculated  for  each  directional  measurement  at  a  specific  location. 

The  second,  or  integral  invariant,  J,  is  given  by 


This  is  a  line  integral  over  the  bounce  path  of  the  particle  and  is  calculated  from  one 
mirror  point  to  the  conjugate  mirror  point  in  the  other  hemisphere  If  one  needs  to 
calculate  the  second  invariant  for  more  than  one  pitch  angle,  more  than  one  integral 
must  be  evaluated.  For  a  given  line  integral  along  a  magnetic  line  of  force  between 
two  conjugate  mirror  points,  many  calls  must  be  made  to  the  magnetic  field 
subroutines.  Since  these  routines  contain  spherical  harmonic  expansions  or  some 
other  equally  complex  expansions  for  the  internal  field,  the  number  of  calls  to  the 
magnetic  field  routines  must  be  minimized.  If  two  particles  have  pitch  angle  a-j  and 
pitch  angle  ap  and  a-|  is  greater  than  ap ,  then  the  bounce  path  length  of  the  pitch 
angle  ap  particle  is  longer  than  that  of  the  «i  particle.  However,  both  particles  will 
follow  the  same  bounce  path  in  the  region  of  overlap.  That  is  the  line  along  which 
the  integral  is  performed  for  pitch  angle  a-]  stops  at  Bmjri  but  that  for  pitch  angle  ap 
continues  on  through  Bmjr-j  to  Bmjr2-  Unfortunately  since  Bmjr  is  inside  the  integral 
sign  the  value  of  the  integral  along  the  line  differs  for  the  two  particles. 
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Maximum  computer  speed  is  obtained  by  dividing  the  invariant  calculation  into  two 
parts.  The  first  part  calculates  the  path  along  the  line  of  force  and  saves  all  of  the 
pertinent  parameters,  and  the  second  part  calculates  the  integral  invariant  for  each 
of  the  pitch  angles.  The  multiple  pitch  angle  invariant  routine,  INVARM,  can 
calculate  the  first  and  second  invariant  of  an  unspecified  number  of  pitch  angles. 

The  angles  must  be  greater  than  zero  (a  pitch  angle  of  zero  would  give  an  infinite 
first  invariant)  and  less  than  or  equal  to  90  degrees.  The  pitch  angle  array  must  be 
sorted  from  biggest  to  smallest  (i.e.  90,80,70,...).  The  line  integral  part  of  INVARM 
steps  along  a  line  of  force  with  a  step  size  that  is  dependent  on  the  curvature  of  the 
line  of  force  until  the  first  Bmax  is  reached.  At  each  step  in  the  integration  the 
program  calculates  the  step  size  as  a  function  of  the  curvature  of  the  field  line.  It 
also  approximates  from  the  present  progress  of  the  integration  the  step  size  needed 
to  reach  Bmax.  It  then  chooses  the  smaller  of  the  two  steps.  It  attempts  to  get 
close  to  Bmax  without  stepping  past  it  on  the  first  approach.  It  is  important  not  to 
exceed  Bmax  since  the  argument  of  the  integral  become  imaginary  if  Bmax  is 
exceeded.  The  step  size  algorithm  appears  to  work  reasonably  well  and  achieves 
an  almost  100%  success  rate  in  not  overstepping  Bmax  on  the  first  try.  If  Bmax  is 
exceeded  the  routine  backs  up  and  attempts  to  determine  a  step  'close'  to  Bmax  but 
smaller  than  Brnax. 

When  the  integration  is  first  started,  the  routine  first  moves  in  the  decreasing  B 
direction  in  order  that  it  can  find  the  precise  value  and  location  of  the  minimum  B. 
When  minimum  B  is  passed,  the  interpolation  routine  determines  a  precise  value  for 
B^in  and  also  determines  the  magnetic  longitude  of  minimum  B.  Once  Bmax  has 
been  found  in  one  direction,  the  integral  is  re-started  at  the  original  location  and  that 
part  of  the  line  to  the  other  mirror  point  is  evaluated.  Once  the  field  line  for  the  first 
pitch  angle  is  found  between  the  two  mirror  points,  the  values  stored  by  the  field  line 

-  9  - 


code  are  used  to  evaluate  the  integral  tor  the  second  invariant  for  the  first  pitch 
angle.  To  calculate  the  invariant  for  additional  pitch  angles,  the  field  line  portion  of 
the  routine  is  reentered  and  the  line  integration  continues  from  the  Bmax  stopping 
point  of  the  previous  pitch  angle.  The  integration  continues  until  the  field  line  up  to 
but  not  exceeding  Bmax  of  the  next  angle  is  determined.  The  integral  for  the 
second  invariant  is  then  calculated.  This  continues  until  the  invariant  of  all  of  the 
pitch  angles  are  determined  or  until  one  of  the  mirror  points,  either  north  or  south  is 
below  1 .03  Re  or  the  maximum  number  of  steps  is  exceeded,  or  until  1 3  Re  is 
exceeded.  Each  subsequent  calculation  utilizes  all  of  the  calculated  values  of  the 
field  strengths  and  step  locations  of  the  previous  pitch  angle,  and  thus  the  number 
of  calls  to  the  magnetic  field  line  routines  is  minimized.  For  example,  the  computer 

time  required  to  calculate  the  invariant  for  18  pitch  angles  (90,  85,  80,  75 . 5)  is 

approximate  2  to  3  times  as  long  as  the  time  required  to  calculate  the  invariant  for 
the  single  pitch  angle  of  20  degrees. 

The  integration  uses  Gill’s  method  of  Runge-Kutta  integration.  This  is  a  fourth  order 
procedure  and  the  error  goes  as  step  size  to  the  fourth  order.  An  internal  error 
control  parameter  can  be  adjusted  to  control  the  errors.  This  parameter  is  set  to 
give  the  "L"  parameter  an  accuracy  of  approximately  0.002.  The  maximum  number 
of  steps  is  100.  Since  it  is  a  fourth  order  procedure,  up  to  400  calls  to  the  magnetic 
field  routines  are  possible.  Typically  on  the  order  10-15  steps  are  required  for  a 
single  pitch  angle  that  mirror  far  off  the  equator  (i.e.  a  pitch  angle  of  20  degrees  at 
an  L  of  5.0).  When  the  invariant  for  18  pitch  angles  are  calculated  an  additional 
step  may  be  needed  at  both  the  northern  and  southern  conjugate  points  for  each 
additional  pitch  angle.  If  successive  pitch  angles  are  very  close  together,  the 
interpolation  routine  may  be  able  to  calculate  the  next  invariant  without  the  need  of 
an  additional  step. 
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When  the  invariant  routine  calculates  the  second  invariant  it  also  integrates  the  total 
column  density  of  the  atmosphere  between  the  mirror  points.  The  density  integral 
uses  the  atmospheric  density  function  developed  for  the  Air  Force  Office  of 
Scientific  Research.  This  function  is  given  by 

density  =  2.7  x  10'1 1  exp[(120-z)/(CON*sqrt(z-103))]  (3) 

Where  z  is  the  altitude  in  kilometer  and  CON  is  an  FI  0.7  dependent  parameter  (70  - 
240)  and  is  given  by 

CON=  0.99  +  0.51 8*  sqrt(F1 0.7/55)  (4) 

Outside  of  3.0  Re  or  outside  of  a  specified  distance,  the  density  function  is  arbitrarily 
set  to  zero,  since  the  function  has  little  validity  above  1000  km  altitude.  It  is, 
however,  a  smoothly  decreasing  function  and  can  thus  be  used  as  an  organizing 
parameter  for  atmospheric  mirror  depth  up  to  3.0  Re. 

Figures  2  -  5  are  copies  of  the  printouts  for  the  calculation  of  the  first  and  second 
invariant,  and  the  'L'  parameter,  and  the  density  totals  for  a  set  of  test  conditions. 
Each  page  has  two  conditions.  The  top  run  uses  the  internal  field  only  and  the 
bottom  run  uses  internal  plus  external  field.  All  runs  are  started  at  latitude  =  0, 
longitude  =1 .0  ,  Day  of  year  =  1  and  Universal  time  =  0.  Figure  2  is  started  at  an 
altitude  of  1 .5  Re.  Both  the  internal  and  internal  plus  external  runs  give  the  same 
result  since  the  external  field  is  unimportant  in  this  region  of  space.  The  small 
variations  in  L  between  the  various  pitch  angle  are  due  to  the  inaccuracies  in  the 
integration  and  more  importantly  to  the  accuracy  and  inherent  approximate 
definition  of  the  L  expansion  (see  Hilton,  J.  Geophys.  Res.  72,  6952,  1971).  Pitch 
angles  smaller  than  35  degrees  have  their  mirror  point  below  200  km  and  the 
invariant  is  thus  not  be  evaluated.  The  internal  field  expansion  diverges  for 
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distances  less  than  1 .0  Re,  and  thus  mirror  points  below  1 .0  cannot  be  assigned  a 
second  invariant.  A  -1 .0  in  any  of  the  parameters  indicates  that  the  mirror  point  is 
too  low  in  the  atmosphere  to  calculate  the  invariant  and  assume  that  the  particle  is 
not  trapped.  A  value  of  100  indicates  an  open  field  line.  As  the  altitude  increases  in 
Figures  3,  4  and  5,  differences  between  the  top  run  with  internal  field  only  and  the 
internal  plus  external  field  run  become  increasingly  large.  The  bottom  run  in  Figure 
5,  a  run  that  calculates  the  invariant  of  a  measurement  at  7.5Re,  shows  that  shell 
splitting  at  7.5  Reis  almost  a  full  Re- 

The  second  to  the  last  column  in  the  printout  shows  the  atmospheric  density 
parameter.  As  the  integration  proceeds  down  the  field  line,  the  integration  sums  the 
atmospheric  density  producing  a  column  number  density  for  the  amount  of 
atmosphere  a  particle  encounters  as  it  bounces  between  the  two  mirror  points.  This 
number  given  in  grams/cm2  is  intended  to  represent  the  importance  of  the 
atmosphere  as  a  loss  mechanism  for  the  trapped  particles  One  observation  can 
already  be  made  from  figure  2;  5  degree  bins  near  the  loss  cone  may  not  be  an 
adequately  fine  resolution.  As  can  be  seen  from  figure  2,  a  5  degree  change  in 
pitch  angle  (40  degrees  to  35  degrees)  causes  the  atmospheric  column  density  to 
change  by  over  two  orders  of  magnitude.  The  next  5  degree  bin  mirrors  below  200 
km.  Figure  3  gives  a  7  order  of  magnitude  change  in  density  encountered  in  the  last 
5  degree  pitch  angle  bin.  These  numbers  are  an  indicator  on  the  expected 
sharpness  of  the  atmospheric  loss  cone 


The  last  column  displayed  in  the  printouts  shows  the  equatorial  pitch  for  the 
specified  particle  at  the  present  location.  It  is  given  by 


aeq  =  sin 


l 


B, 


-sin(«llvjl) 


(5) 
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where  aeq  is  the  equatorial  pitch  angle  and  a|0ca|  is  the  local  pitch  angle.  It  is 
important  to  remember  that  the  equatorial  pitch  angle  of  a  particle  is  not  an 
invariant.  The  equatorial  pitch  angle  of  a  particle  changes  as  the  particle  drifts 
around  the  earth.  This  effect  is  most  important  at  large  distances  where  the  earth’s 
field  becomes  more  asymmetric. 

INVARM,  efficiently  calculates  the  first  and  second  invariant  of  numerous  pitch 
angles  at  a  single  satellite  position.  It  also  calculates  the  effective  L’s  for  the  given 
set  of  invariant.  It,  furthermore,  determines  the  actual  minimum  B  field  along  the 
field  line,  the  magnetic  latitude  of  the  observation  point  and  the  magnetic  longitude 
where  the  field  line  crosses  the  magnetic  equator.  It  also  provides  the  column 
density  of  the  atmosphere  between  the  mirror  points  and  thus  is  able  to  estimate  the 
amount  of  scattering  or  absorption  that  can  take  place  at  the  specified  pitch  angle. 
INVARM  provides  a  complete  characterization  of  all  of  the  pertinent  magnetic 
parameters  for  any  set  of  pitch  angles  anywhere  within  the  stable  trapping  region. 

Appendix  C  lists  routine  INVARM,  a  test  routines,  and  the  various  subroutines  and 
functions  required  for  its  correct  operation.  The  magnetic  field  routines  required  to 
define  the  field  are  described  in  Appendix  A  and  B. 

3.0  Summary 

The  new  pitch  angle  dependent  invariant  code  with  its  ability  to  calculate 
atmospheric  mirror  point  depth  will  substantially  improve  our  ability  to  organize  the 
directional  charged  particle  data  from  the  CRRES  satellite  and  develop  the  next 
generation  radiation  belt  models.  Considerable  effort  was  expended  to  achieve  an 
efficient  code  that  minimizes  computer  time.  The  speed  of  the  IGRF  internal  model 
was  considerably  improved  and  the  invariant  code  is  highly  optimized.  The  success 
of  this  code  must  now  be  verified  by  organizing  the  CRRES  data  with  this  code. 
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Figure  2 
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Figure  3 
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Figure  4 
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Appendix  A 

Internal  Magnetic  Field  Subroutines 

Subroutine  SPIGRF  is  a  fast  version  of  the  IGRF  internal  field  subroutine.  Instead  of  using 
DO  loops  to  expand  the  spherical  harmonic  coefficients,  it  writes  out  the  expansion 
according  to  a  routine  developed  by  J.  Cain.  This  version  of  the  routine  developed  for  the 
CRRES  program  is  not  designed  for  stand  alone  use  and  is  designed  to  be  a  part  of  a  total 
magnetic  field  program  that  includes  the  internal  as  well  as  the  external  magnetic  field. 

The  calling  arguments  are  thus  passed  in  a  common  block  and  the  geomagnetic  latitude 
and  longitude  are  passed  via  their  sines  and  cosines.  This  techniques  saves  computer 
time.  A  stand-alone  version  would  be  much  easier  to  use  if  the  calling  variables  were 
transmitted  as  standard  subroutine  arguments.  If  an  easy  to  use  stand  alone  version  is 
needed,  a  simple  change  to  the  first  few  lines  of  the  code  can  produce  a  very  efficient 
stand-alone  internal  field  code. 

The  method  of  coding  the  magnetic  field  in  SPIGRF  is  inherently  faster  than  a  DO  loop 
version.  Furthermore,  the  addition  of  a  term  dropping  algorithm  increases  the  speed  as 
much  as  a  factor  of  35.  The  variable  CONA  dimensioned  1 1  contains  the  altitudes  at 
which  successive  terms  are  to  be  dropped.  A  linear  interpolation  between  these  distance 
values  drops  the  terms  off  smoothly.  The  smooth  feature  is  important  since  it  prevents 
discontinuities  in  the  magnetic  field  from  disrupting  the  integration  steps  and  the 
interpolation  algorithms  in  the  field  line  tracing  program. 

A.1  Calling  Sequence 

The  transfer  of  information  between  SPIGRF  and  the  calling  routine  is  performed  via 
labeled  COMMON  GCOM. 

INPUT  values 

YEARI  Contains  the  year  for  which  the  coefficients  are  to  be  determined.  The 

supplied  coefficients  are  valid  from  1945  to  the  present.  The  1945  coefficients 
are  used  for  years  earlier  than  1 945.  Predicting  far  into  the  future  is 
hazardous  since  the  time  derivative  terms  do  not  have  long  term  validity.  If 
YEARI  is  changed  by  .1  years  since  the  last  call  a  new  updated  set  of 
coefficients  is  calculated.  It  is  suggested  that  YEARI  be  used  to  set  the 
desired  epoch  and  then  left  constant. 

NMAXN  Contains  the  number  of  terms  desired.  If  this  is  left  0  then  the  full  1 1  term 
expansion  of  IGRF  is  used.  If  NMAXN  is  between  2  and  10,  then  the 
maximum  number  of  terms  is  set  to  that  number.  Term  dropping  still  takes 
place  for  larger  distances. 

ST  Sine  of  the  geographic  co-latitude 

CT  Cosine  of  the  geographic  co-latitude 
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SPH  Sine  of  the  geographic  longitude. 

CPH  Cosine  of  the  geographic  longitude 

AOR  6371 ,2/R,  where  R  is  the  distance  from  the  center  of  the  Earth  in  km 

OUTPUT  values 

BR  The  radial  component  of  the  magnetic  field  in  gauss  (ie.  nanotesla) 

BT  Theta  component  (south  pointing)  component 

BP  Phi  component  (East) 

Each  time  the  field  coefficients  are  updated,  the  value  of  the  new  dipole  moment  is  stored 
in  labeled  COMMON  /MOMENT/  XM.  It  is  thus  available  for  use  by  any 
routine  that  needs  it,  such  as  the  L  value  program. 


NOTE:  SPIGRF  uses  a  true  spherical  coordinate  system  with  the  z  axis  along  the 

geographic  north  pole,  the  x  axis  through  the  longitude  of  Greenwich.  R, 
Theta  and  Phi  are  the  true  spherical  polar  coordinates. 


The  routine  will  read  several  of  the  IGRF  Coefficient  sets.  The  coefficient  sets  are  listed  at 
the  end  of  this  appendix.  Subroutine  FLDCOF  sets  the  FORTRAN  logical  unit  for  reading 
the  coefficients  to  1 1 .  The  actual  read  statement  for  reading  the  coefficients  are  found  in 
subroutine  GETGAU. 
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SUBROUTINE  SPIGRF 


VERSION  4/91 

WRITTEN  BY  K.  A.  PFITZER  (714)  896-3231 

SPIGRF  IS  A  MODIFIED  VERSION  OF  J.C.  CAIN'S  14  TERM  FAST  SPHRC 
ROUTINE. 

IT  HAS  BEEN  SHORTENED  TO  11  TERMS  FOR  CONSISTENCY  WITH  THE  IGRF 
COEFFICIENT  SET. 

IT  HAS  A  TRUNCATION  FOR  LARGE  R  -  THE  TRUNCATION  BETWEEN  TERMS 
IS  SMOOTH  AND  MAINTAINS  AN  ACCURACY  OVER  THE  NON-TRUNCATED  VERSION 
OF  BETTER  THAN  .1  NANOTESLA. 

DEPENDING  ON  ALTITUDE  THIS  VERSION  RUNS  FROM  1.5  TO  35.0  TIMES  AS  FAST 
AS  THE  STANDARD  SCHMITT  NORMLIZED  IGRF  ROUTINES 
THE  SUPPORT  ROUTINES  READ  THE  STANDARD  IGRF  COEFFICIENTS  AND 
CONVERT  THEM  TO  GAUSS  NORMALIZED  FOR  USE  BY  THIS  ROUTINE 


The  first  time  the  routine  is  called,  the  routine  calls  routine 
call  routine  FLDCOF  to  obtain  the  correct  IGRF  coefficients.  If 
the  date  changes  by  more  than  .1  year  the  coefficients  are  updated, 
new  coefficients  are  obtained  if  required. 


INPUT  -- 
YEAR  I 

ST 

CT 

SPH 

CPH 

AOR 

NMAXN 


OUTPUT  - 
BR 
BT 
BP 


COMMON  BLOCK  GCOM 

IS  THE  YEAR,  IF  YEARI  CHANGES,  THE  COEFFICIENTS  ARE 
UPDATED . 

SINE  OF  THE  GEOGRAPHIC  CO-LATITUDE. 

COSINE  OF  THE  GEOGRAPHIC  CO-LATITUDE. 

SINE  OF  THE  GEOGRAPHIC  LONGITUDE. 

COSINE  OF  THE  GEOGRAPHIC  LONGITUDE. 

6371. 2/R,  WHERE  R  IS  THE  GEOCENTRIC  DISTANCE  IN  KM  FROM 
THE  CENTER  OF  THE  EARTH. 

MAXIMUM  NUMBER  OF  TERMS  TO  BE  USED  (MUST  BE  LESS  OR 
EQUAL  TO  11) .  THIS  ROUTINE  PRESETS  IT  TO  11 
NMAXN  OF  11  CORRESPONDS  TO  THE  10TH  ORDER  IGRF  MODELS 
IF  NMAXN  IS  >2  AND  <11,  NMAXN  TERMS  ARE  USED,  ELSE  THE 
NUMBER  OF  TERMS  USED  IS  11  OR  THE  MAXIMUM  TERMS  IN  THE 
IGRF  DATA  SET. 

-  COMMON  BLOCK  GCOM 

RADIAL  COMPONENT  OF  FIELD  IN  GAUSS. 

THETA  COMPONENT  (SOUTH  POINTING)  COMPONENT. 

PHI  COMPONENT  (EAST) 


DIMENSION  G(ll, 11) , CONST ( 1 1 , 1 1 ) , FM ( 1 1 ) , FN ( 1 1 ) 

COMMON  /MODEL/G 

COMMON  /GCOM/  ST , CT , SPH, CPH, AOR, BT, BP , BR, NMAXN, YEARI 
COMMON  /MOMENT /XM 
DIMENSION  CONA (11) 

DATA  YRLAST  /-12345./ 

DATA  IFIRST/0/ 

DATA  CONA/ 0 . , 12.0,8.0,6.0,5.0,4.0,3.2,2.5,2.0,1.6,1.4/ 


SET  UP  INITIAL  CONSTANTS  DURING  FIRST  CALL 

IF { IFIRST . NE . 0)  GO  TO  199 

IFIRST=1 

FM ( 1 ) =0 

DO  6  N=2 , 1 1 

FM (N) =N- 1 

FN (N) =N 

DO  6  M= 1 , N 

6  CONST (N,M) =FLOAT ( (N-2)  **2-(M~i) **2)/FLOAT( ( 2  *N- 3 )  *  (2*N-5)  ) 


SET  UP  THE  COEFFICIENTS 

IF  YEARI  HAS  CHANGED  BY  MORE  THAN  .1  YEAR  UPDATE  THE  COEFFICIENTS 
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199 


IF (ABS (YRLAST-YEARI) . LT . 0 . 1 )  GO  TO  230 
CALL  FLDCOF ( YEARI , DIMO, MAXN) 

XM=D IMO/ 1 . 0E5 
YRLAST=YEAR I 


230  NMAX=MAXN 

I F ( NMAXN . GE . 2 . AND . NMAXN . LT . MAXN ) NMAX=NMAXN 

AR=AOR*AOR*AOR 

C2=G (2,2) *CPH+G (1,2) *SPH 

ER=- (AR+AR) * (G (2, 1) *CT+C2*ST) 

BT=AR* (C2*CT-G (2, 1) *ST) 

BP=AR* (G { 1 , 2 ) *CPH-G (2,2) *SPH) 

IF  (NMAX.LE.2)  RETURN 
R=1 . /AOR 

IF (R . GT . CONA ( 2 ) )  RETURN 
CON=0 . 

SP2=SPH 

CP2  =  CPH 

P2 1 =CT 

P22=ST 

DP2l=-ST 

DP22=CT 

N=3 


1-00006 

1-00007 

1-00008 

1-00009 

1-00010 

1-00011 


1-00012 

1-00013 

1-00014 

1-00015 

1-00016 

1-00017 


SP3=  (SP2+SP2) *CP2 

CP 3= (CP2  +  SP2) *  (CP2-SP2 ) 

P31=CT*P2i-CONST (3, 1) 

P32=CT*P22 

P33=ST*P22 

DP31=-P32-P32 

DP32=CT*DP22-P33 

DP33=-DP31 

C2=G (3, 2) *CP2+G (1, 3) *SP2 
C3=G (3,3) *CP3+G (2,3) *SP3 
AR=AOR*AR 

XR=BR-FN  (3)  *  AR*  (G  (3,  1)  *P31+C2*P32+C3*P33) 

XT=BT+AR* (G ( 3, 1) *DP31+C2  *DP32+C3*DP33) 

XP=BP - AR* (FM(2) * (G(3, 2) *SP2-G(1, 3) *CP2) *P21+FM(3) * 
+3) *CP3) *P22 ) 

BP=BP  *ST 
XP=XP  *ST 

IF (NMAX . LE . 3 )  GO  TO  21 
IF ( R . GT . CONA ( 3 ) )  GO  TO  20 
N=  4 


1-00019 

1-00020 

1-00021 

1-00022 

1-00023 

1-00024 

1-00025 

1-00026 

1-00027 

1-00028 

1-00029 

1-00030 

1-00031 

(G(3,  3)  *SP3-G(2, 1-00032 
1-00033 
1-00035 
1-00035 


SP4=SPH*CP3+CPH*SP3 

CP4=CPH*CP3-5PH*SP3 

P41=CT*P3 1 -CONST ( 4 , 1} *P21 

DP41=CT*DP31-ST*P31 -CONST ( 4 , 1) *DP21 

P42=CT*P 32 -CONST  (4,2)  *F22 

DP42=CT*DP  32 -ST *P 32 -CONST (4,2) *DP22 

P43=CT*P33 

DP43=CT*DP33-ST*P33 

P44=ST*P33 

DP44  =  FM (4)  *P4  3 

C2=G ( 4  ,  2)  *CP2+G(1, 4)  *SP2 

C3=G (4, 3) *CP3+G (2, 4) *SP3 

C4=G(4, 4) *CP4+G(3, 4) *SP4 

AR=AOR*AR 

BR=XR-FN ( 4 ) * AR* (G ( 4 , 1 ) * P 4 1 +C2 *P42+C3 *P4 3+C4 *P4 4 ) 
BT=XT+AR* ( G ( 4 , 1) *DP41+C2*DP42+C3*DP43+C4*DP44) 
BP=XP-AR* (FM (2) * (G (4, 2) »SP2-G ( 1, 4) *CP2) *P42+FM (3) 
M)  *CP 3 )  *P43*-FM(4)  *  (G(4,  4)  *SP4-G(3,  4)  *CP4)  *P4  4) 

IF (NMAX. LE . 4)  GC  TO  11 
IF ( R . GT . CON A ( 4 ) )  GO  TO  1 0 
N=  5 


1-00037 

1-00038 

1-00039 

1-00040 

1-00041 

1-00042 

1-00043 

1-00044 

1-00045 

1-00046 

1-00047 

1-00048 

1-00049 

1-00050 

1-00051 

1-00052 

(G (4, 3) *SP3-G (2, 1-00053 
1-00054 
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SP5= (SP3+SP3) *CP3 
CP5= (CP3+SP3) * (CP3-SP3) 

P  51=CT*P4 1 -CONST ( 5 ,  1) *P31 

DP5 1=CT*DP4 1 -ST*P  4 1-CONST ( 5 , 1) *DP31 

P52=CT*P42-CONST<5,2)  *P32 

DP52=CT*DP42 -ST*P42-CONST (5,2) *DP32 

P53=CT*P43-CONST (5, 3) *P33 

DP53=CT*DP43-ST*P43-CONST (5, 3) *DP33 

P54=CT*P44 

DP54=CT*DP44-ST*P44 

P55=ST*P44 

DP55=FM ( 5) *P54 

C2=G (5,2) *CP2+G (1,5) *SP2 

C3=G ( 5, 3) *CP3+G (2,5) MP3 

C4=G (5,4) *CP4+G (3, 5)  MP4 

C5=G (5, 5) *CP5+G (4,5) MP5 

AR=AOR*AR 

XR=BR-FN(5) *AR* (G(S, 1) *P5 1 +C2 *P 52 +C3 *P53+C4 *P5 4 +C5 * P 55 ) 

XT=BT+AR* (G (5, 1) *DP5 1 +C2 *DP52+C3*DP 53+C4 *DP 5 4+C5 *DP 5 5 ) 

XP=BP-AR*  (FM  (2 )  *  (G  ( 5 , 2 )  MP2-G  (1, 5)  *CP2)  *P52+FM(3)  MG  (5,  3)  MP3-G 
+  5)  *CP3)  *P53+FM  (  4 )  *  (G(5,  4)  MP4-G(3,  5)  *CP4) *P54+FM (5)  MG(5,  5)  MP5- 
+4,5) *CP5 ) *P55) 

IF (NMAX . LE . 5 )  GO  TO  21 
IF (R.GT .CONA (5) )  GO  TO  20 
N=6 

SP6=SPH*CP5+CPHMP5 

CP6=CPH*CP5-SPHMP5 

P61=CT*P51 -CONST (6, 1) *P41 

DP 61=CT*DP5 1-ST*P51 -CONST (6, 1) *DP41 

P62=CT*P 52 -CONST (6,2) *P42 

DP62=CT*DP52-ST*P52-CONST (6, 2) *DP42 

P63=CT*P53-CONST (6,3) *P43 

DP63=CT*DP53-ST*P53-CONST (6, 3) *DP43 

P64=CT*P5 4-CONST (6,4) *P44 

DP 64=CT*DP54 -ST *P 54-CONST (6, 4) *DP44 

P65=CT*P55 

DP6S=CT*DP55-ST*P55 

P66=ST*P55 

DPb6=FM(6) *P65 

C2=G (6,2) *CP2+G ( 1 , 6) *SP2 

C3=G  (6,3)  *CP3+G  (2  ,  6)  MP3 

C4=G (6,4) *CP4+G (3,6)  *SP4 

C5=G ( 6, 5) *CP5+G(4, 6) *SP5 

C6=G ( 6, 6) *CP6+G { 5, 6) ‘S26 

AR=AOR*AR 

BR=XR-FN (6) *AR* (G (6, 1) * P 61 +C2 *P 62+C3 *P 63+C4 *P 64 rC5 *P 65+C 6  * P 66 ) 
BT=XT  +AR* (G ( 6, 1) *DP 6 1 +C2 *DP 62+C3 *DP 63+C4 *DP 64 +C5 *DP 65+C6 *DP 66 ) 
BP=XP-AR*  (FM(2)  *  (G  ( 6, 2)  *SP2-G(1, 6)  *CP2)  *P62+FM(3)  MG  (  6,  3)  *SP3-G  ( 
+  6)  *CP3)  *P  63  +  FM  (  4  )  *  (G  (  6,  4 )  *SP  4  -G  (3,  6)  *CP4)  *P64+FM(5)  MG(6,5)'SP5- 
+4,6) *CP5) *P65+FM ( 6) » (G(6, 6) *SP6-G(5, 6) *CP6) *P66) 

IF (NMAX. LE. 6)  GO  TO  11 
IF ( R . GT . CONA ( 6 ) )  GO  TO  10 
N=7 

SP7=  (SP4  +  SP4 )  *CP  4 

CP 7=  (CP4+SP4 )  *  (CP4-SP4) 

P71=CT*P6 1 -CONST (7, 1 ) *P  5 1 
DP71=CT*DP61-ST*P61 -CONST (7,1) *DP51 
P72=CT*P 62 -CONST (7, 2) *P52 
DP72=CT*DP62-ST*P 62 -CONST (7,2) *DP52 
P73=CT*P 63 -CONST (7, 3) *P53 
DP73=CT*DP63-ST*P 63-CONST (7, 3) *DP53 
P74=CT*P64-CONST (7, 4) *°54 
DP74=CT*DP64-ST*P 64 -CONST (7, 4) *DP54 
P75=CT*P 65 -CONST (7, 5) *P55 


-00057 

-00058 

-00059 

-00060 

-00C61 

-00062 

-00063 

-00064 

-00065 

-00066 
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-00068 

-00069 

-00070 

-00071 

-00072 
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-00078 


-00081 

-00082 

-00083 

-00084 

-00085 

-00086 

-00087 

-00088 

-00089 

-00090 
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-00092 
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-00094 
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-00099 

-00100 
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-00102 

-00103 

-00104 

-00105 


-ooice 

-00109 

-00110 

-00111 

-00112 

-00113 

-00114 

-00115 

-00116 

-00117 

-00118 
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DP7  5=CT  *DP  65-ST  *P  65 -CONST (7,5) *DP55 

1-00119 

P7  6=CT*P66 

1-00120 

DP 7  6=CT*DP  66-ST*P66 

1-00121 

P77=ST*P66 

1-00122 

DP7  7  =  FM ( 7 ) *P7  6 

1-00123 

C2=G (7,2) *CP2+G (1,7) *SP2 

1-00124 

C3=G  (7,3)  *CP3  *-G  (2 , 7  )  *SP3 

1-00125 

C4=G (7, 4) *CP4+G (3, 7) *SP4 

1-00126 

C5  =  G (7,5) *CP  5+G ( 4 , 7 )  *  SP  5 

1-00127 

C  6=G ( 7 , 6) *CP  6+G  (5,7)  *SP6 

1-00128 

C7=G (7, 7) *CP7+G (6, 7) *SP7 

1-00129 

AR=AOR*AR 

1-00130 

XR=BR-FN ( 7 ) *AR* (G(7, 1) *P7 1 +C2 *P72 +C3 *P7 3+C4 

*P7  4+C5  *P7  5+C6*P7  6+C7  *P 1-0013 1 

+  77) 

1-00132 

XT  =  BT  +AR* (G (7, 1)  *DP7 1 +C2 *DP72+C3 *DP7 3+C4 *DP7 4 +C5 *DP75+C6*DP7 6+C 7 *D1 -001 33 
+P77 )  1-00134 

XP=BP-AR* (FM (2) *(G<7,2) *SP2-G(1, 7) *CP2) *P72+FM<3) * (G (7, 3) *SP3-G(2, 1-00135 
+  7 ) *CP  3 ) *P73  +  FM (4) * (G(7, 4) *SP4-G(3,7) *CP4) *P74+FM(5) * (G ( 7 , 5) *SP5-G ( 1-00136 
+  4, 7)  *CP5) *P7  5  +  FM (6) * (G(7, 6) *SP6-G(5,7) *CP6) *P76+FM(7) * (G (7, 7) *SP7 -1-001 37 
+G ( 6, 7) *CP7 ) *P77 )  1-00138 

IF (NMAX. LE . 7)  GO  TO  21 
IF ( R . GT . CON A ( 7 ) )  GO  TO  20 
N  =  8 


SP8=SPH*CP7+CPH*SP7 
CP8=CPH*CP7-SPH*SP7 
P81=CT*P7 1 -CONST (8, 1) *P61 
DPS 1=CT*DP71-ST*P7 1 -CONST (8,  1) *DP61 
P82=CT*P 72 -CONST (8, 2)  *P62 
DP82=CT*DP72-ST*P72-CONST  <8, 2) *DP62 
P83=CT*P73-CONST (8, 3) *P63 
DP83=CT*DP73-ST*P7 3-CONST (8, 3) *DP63 
P84=CT*P74-CONST (8,4) *P64 
DP8  4=CT  *DP7  4 -ST  *P7  4-CONST ( 8 , 4)  *DP64 
P85=CT*P75-CONST (8,5) *P65 
DP85=CT*DP75-ST*P75-CONST (8, 5) *DP65 
P86=CT*P7 6-CONST (8,  6) *P66 
DP36=CT*DP76-ST*P7 6-CONST (8, 6) *DP66 
P87=CT  *P7  7 
DP87=CT*DP77-ST*P77 
P33=ST*P77 
DP88=FM (8) *P87 
C2=G ( 8 , 2 ) *CP2  +G ( 1 ,  8 ) 

C3=G (8, 3) *CP3+G (2,8) 

C4=G (8, 4)  *CF  4+G (3, 8) 

C5=G ( 8 , 5) *  CP5+G ( 4 , 8) 

C6=G  (8, 6) *CP6+G (5, 8) 

C7=G  (8,7)  *CP7 +G ( 6, 8) 

C  8  =  G ( 9 , 8)  *CF 8-0(7, 8) 

AR=AOR  *  AR 


-SP2 
‘SP3 
'SP4 
'  SP5 
1  SP6 
'  SP7 
‘SF  8 


BR=XR-FN ( 8 ) * AR* (G(8, 1) * P8 1 +C2 *P 82+C3 ‘P83+C4 *P 84+C5 *P85+C6*P8 6+C7 *P 1 

*  87  +C8  *P88 )  1 

BT  =  XT+AR  * (G (8,1) *I'P81  +  C2  *DP82 tC3‘DP83+C4*DP84+C5‘DP85+C6*DP86+C7*Dl 

*  P67  +C9*DP«-:  )  1 

FiP -XP-AR*  (FM  !  2!  *  (G(«,2)  *St  2-G  (1 , 8)  *CP2)  *P82+FM  (3)  *(G(8,3)*SP3-G(2,1 

*  9)  *CP3)  ‘P83-FM  (4)  *  (G  (8,  4)  *SP4-G  (3,  8)  *CP4)  *P8  4+FM  (5)  *  (G(8,  5)  *SP5-G  (1 

•4,8)  *CP 5 )  *F  8  5  *■  FM  (  6 )  *  (G  ( 8 ,  6)  *SP6-G(5,  8)  *CP6)  *P86+FM(7)  *  (G(8,  7)  *SP7-1 
•G  <6, 8)  *CP7)  *P87  *FM (6)  *  <G (8, 8)  *SP8-G(7, 8)  *CP8) *P88)  1' 

IF (NMAX. LE. 9)  GO  TO  11 
IF (R .GT .CGNA (8) )  GO  TO  10 
N  -  9 

SP9- (SP5+SP5) *CP5  1 

CP (CP 5  +  C.r  1  )  *  (CP5-SP5)  1 

F  9 1 =CT *  P  8 1 - CC NOT ( 0, 1 )  *  P7 1  1 

DF 0 1=CT  *OP  9  )  -''T*P 31  -CONST  (9,1)  *DP7 1  1 


-00141 

-00142 

-00143 

-00144 

-00145 

-00146 

-00147 

-00148 

-00149 

-00150 

-00151 

-00152 

-00153 

-00154 

-00155 

-00156 

-00157 

-00158 

-00159 

-00160 

-00161 

-00162 

-00163 

-00164 

-00165 

-00166 

-00167 

-00168 

-00169 

-00170 

-00171 

-00172 

-00173 

-00174 


•00177 

•00178 

•00179 

•00180 
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P92=CT*P82-CONST <9, 2) *P72  1-00181 

DP92=CT*DP82-ST*P82-CONST (9,2) *DP72  1-00182 

P 9 3=CT*P 83-CONST ( 9, 3) *P7  3  1-00183 

DP93=CT*DP83-ST*P83-CONST (9, 3) *DP73  1-00184 

P94=CT*P84-CONST ( 9, 4) *P74  1-00185 

DP94=CT*DP84-ST*P84-CONST { 9, 4) *DP74  1-00186 

P95=CT*P 85-CONST ( 9, 5) *P75  1-00187 

DP95=CT*DP85-ST*P85-CONST (9, 5) *DP75  1-00188 

P96=CT*P8 6-CONST (9, 6) *P76  1-00189 

DP96=CT*DP86-ST*P86-CONST (9, 6) *DP76  1-00190 

P97=CT*P 87 -CONST (9,7) *P77  1-00191 

DP  97=CT*DP87 -ST*P87 -CONST (9,7) *DP77  1-00192 

P98=CT*P88  1-00193 

DP98=CT*DP88-ST*P88  1-00194 

P99=ST*P88  1-00195 

DP99=FM ( 9) *P98  1-00196 

C2=G ( 9, 2) *CP2+G (1,9) *SP2  1-00197 

C3=G ( 9, 3) *CP3+G (2 , 9) *SP3  1-00198 

C4=G ( 9, 4 ) *CP4+G (3, 9) *SP4  1-00199 

C5=G (9, 5) *CP5+G(4, 9) *SP5  1-00200 

C6=G (9, 6) *CP6+G ( 5, 9) *SP6  1-00201 

C7=G ( 9, 7 ) *CP7+G ( 6, 9) *SP7  1-00202 

C8=G (9, 8) *CP8+G (7, 9) *SP8  1-00203 

C9=G(9, 9) *CP9+G(8, 9) *SP9  1-00204 

AR=AOR*AR  1-00205 

XR=BR-FN ( 9) *AR* (G(9, 1) *P 91+C2 *P 92+C3 *P 93+C4 "P94+C5 *P 95+C6 *P 96+C7 *P 1 -002 0 6 
+97+C8*P98+C9*P99)  1-00207 

XT=BT+AR* (G ( 9, 1 ) *DP91+C2*DP  92+C3*DP93+C4  *DP  94+C5  *DP  95+C  6*DP96+C7*D1-00208 
+P97+C8*DP98+C9*DP99)  1-00209 

XP=BP-AR* (FM (2 ) * (G ( 9, 2 ) *SP2-G<1, 9) *CP2) *P92+FM(3) * (G (9, 3) *SP3-G (2, 1-00210 
+9) *CP3) *P93+FM(4) *<G(9, 4) *SP4-G{3, 9) *CP4) *P94+FM(5) * (G(9, 5) *SP5-G ( 1-00211 
+4, 9) *CP5) *P95+FM ( 6) * (G(9, 6) *SP6-G{5, 9) *CP6) *P96+FM(7) * (G (9, 7) *SP7 -1-002 12 
+G ( 6, 9) *CP7 ) *P97+FM (8) * (G(9, 8) *SP8-G(7, 9) *CP8) *P98+FM(9) * (G (9, 9) *SP 1-002 13 

1-00214 


+9-G (8, 9) *CP9) *P99) 
IF (NMAX . LE . 9 )  GO  TO 
IF (R . GT . CONA ( 9) )  GO 
N=10 


21 

TO 


20 


SP10=SPH*CP9+CPH*SP9 

CP10=CPH*CP9-SPH*SP9 

P 1 0 1=CT*P  91-CONST (10, 1) *P81 

DP101=CT*DP91-ST*P91-CONST (10, 1) *DP81 

P102=CT*P92-CONST  (10,2) *P82 

DP102=CT*DP92-ST*P 92 -CONST (10,2) *DP82 

P103=CT*P 93-CONST (10, 3) *P83 

DP  1 03=CT*DP 93 -ST *P 93-CONST (10, 3) *DP83 

P104=CT*  P  94 -CONST (10, 4 )  *P84 

DP104=CT*DP94-ST*P94-CONST (10, 4) *DP84 

P105=CT*P 95 -CONST (10, 5) *P85 

DP105=CT*DP95-ST*P95-CONST (10, 5) *DP85 

P106=CT*P 96-CONST (10, 6) *P86 

DP 106=CT*DP96-ST*P 96-CONST (10, 6) *DP86 

P 1 07=CT*P  97 -CONST (10,7) *P87 

DP107=CT*DP97-ST*P97-CONST  (10,7) *DP87 

P108=CT*P 98-CONST (10, 8) *P88 

DP 108=CT*DP98-ST*P 98 -CONST (10, 8) *DP88 

P109=CT*P99 

DP  1 0  9=CT*DP  99-5T  *P  99 

P1010=ST*P99 

DP  10 1 0  =  FM (10) *P109 

C2=G (10,2) *CP2+G ( 1 , 10) *SP2 

C3=G (10, 3) *CP3+G (2 , 10) *SP3 

C4=G (10, 4 ) *CP4  +G ( 3 , 10) *SP4 

C5=G ( 10, 5) *CP5+G(4, 10) *SP5 

C6=G (10,6) *CP6+G (5, 10) *SP6 


1-00217 
1-00216 
1-00219 
1-00220 
1-00221 
1-00222 
1-00223 
1-00224 
1-00225 
1-00226 
1-00227 
1-00228 
1-00229 
1-00230 
1-00231 
1-00232 
1-00233 
1-00234 
1-C0235 
1-00236 
1-00237 
1-00238 
1-00239 
1-00240 
1-00241 
1  -00242 
1-00243 


-  24  - 


C7=G (10, 7 ) *CP7+G ( 6, 10) *SP7  1' 

C8=G(10, 8) *CP8fG(7, 10) *SP8  1' 

C9=G ( 1 0, 9) ‘CP9+G ( 8, 10) *SP9  1 

C10=G  ( 10,  10)  *CP10+-G(9,  10)  *SP10  1 

AR=AOR*AR  1 

BR=XR-FN  ( 10  )  *  AR  *  (GUO,  1)  *P  1 01+C2  *P  102  +C3‘P  103+C4  *P  1 04  +C5  *P1 05+C6*P  1 
+106+C7*P 107+C8 *P108+C9*P1094C10*P10 10)  1 

BT=XT+AR* (G (10, 1 ) *DP101+C2*DP102+C3*DP103+C4*DP104+C5*DP105+C6*DP11 
‘06+C7*DP107+C8*DP108+C9*DP109+C10  * DP  10 1 0 )  1 

BP=XP-AR*  (FM  (2)  MG  (10,  2)  *SP2-G  (1,  10)  *CP2)  *P102+FM(3)  *  (G  (10,  3)  *SP3-1 
+G  (2, 10) *CP3)  *P 1 03  +FM ( 4 ) * (G ( 10, 4 ) *SP4-G (3, 10) *CP4) ‘P104+FM (5) MG  (101 
+, 5) ‘SP5-G (4, 10) *CP5) ‘P105+FM (6) MG (10, 6) *SP6-G(5, 10) *CP6) *P106+FM(1 
+  7) *  (G ( 10 , 7 ) *SP7-G(6, 10) *CP7 ) *P107+FM(8) *(G(10,8)*SP8-G(7,10) *CP8) *1 
+P108+FM  ( 9)  *  (GUO,  9)  *SP  9-G ( 8 ,  10)  *CP9)  *P109+FM(10)  MG  ( 10,  10)  *SP10-G  ( 1 
+9, 10) *CP10) *P1010)  1 

IF(NMAX.LE.IO)  GO  TO  11 
IF (R.GT.CONA(IO) )  GO  TO  10 

N=1 1  1' 

SPU=  (SP6  +  SP6)  *CP6  1 

CPU=  (CP6  +  SP6)  *  (CP6-SP6)  1 

P 1 1 1=CT*P1 0 1 -CONST (11, 1) *P91  1 

DP111=CT*DP101-ST*P101 -CONST (11, 1) *DP91  1 

P 1 12=CT*P 102-CONST (11,2) *P92  1 

DPI 12=CT  *DP 1 02 -ST  *P 1 02 -CONST (11,2) *DP92  1 

P 1 1  3=CT*P 1 03-CONST ( 11 , 3) *P93  1 

DP 113=CT*DP103-ST*P 103-CONST (11, 3) *DP93  1 

PI  1 4=CT*P104 -CONST (11, 4 ) *P94  1 

DPI 1 4=CT  *DP 104-ST*P 104 -CONST (11, 4 ) *DP94  1 

P 1 15=CT*P 105-CONST (11, 5) *P95  1 

DP  1 1  5=CT*DP  1 05-ST ‘P105-CONST (11,5) *DP95  1 

PI 1 6 =CT*P 10  6-CONST (11, 6) *P96  1 

DP  116=CT*DP106-ST*P 10 6-CONST (11,  6)  *DP96  1 

P 1 1 7=CT  *P 1 07 -CONST (11,7)  *P  97  1 

DP117=CT*DP107-ST*P 1 07 -CONST (11,7) *DP97  1 

P 1 1 8 =CT*P1 08-CONST (11, 8 ) *P98  1 

DPI 18=CT*DP108-ST*P 108-CONST (11, 8) *DP98  1 

P119=CT*P109-CONST (11, 9) *P99  1 

DP U9=CT*DP10  9-ST*P  10 9-CONST  (11, 9)  *DP99  1 

P 1 1 10=CT*P 1 010  1' 

DP  1 1 1 0=CT*DP 1010-ST*P1010  1' 

P  1 1 1  =  ST*P  1  0 1  0  1' 

DPI 1 1 1=FM (II) *P1 1 10  1- 

C2=G(11,2) ‘CP2+G(1, 11) *SP2  1 

C3=G ( 1 1 , 3) *CP3+G(2, 11) *SP3  l- 

C4=G ( 1 1 , 4 )  *C?4+G ( 3 ,  11) *  SP4  1- 

C5=G (11, 5)  *  C  P  5  +  G ( 4 , 11) *SP5  1' 

C6-G ( 1 1 , 6) *  CP  6  4  G ( 5, 11) *SP6  1 

C7=G (11,7)  *  P7 +0(6, 11) *SP7  1 

C8-GU1, 8)  *  C P  8  +  0(7,  ]  1)  *SP8  1 

C9=G (11, 9>  *CP9+G (8, 11) *SP9  1- 

Cl 0=G (11, 10) *CP10»G(9, 11) *SP10  1- 

ci  i  =a (ii,ii)  *cm  +guo,  id  *spii  !■ 

AF<=-ACR*AR  1 

BP  =  BK-FN ( 1  1)  *  AR  * (G(ll,  1) *F  1  1 1 +C2 *P 1 12+C3 * P 1 1 3+C4 * P 1 1 4 +C5 *P 1 1 5+C6*P 1 
U16+C7*P1 1  7+C8*Pl  i8-tC9*P119+C10*P1110+CU*Pllll)  1- 

BT-BT  +  AR* (G (11,1)  'DPI  1 1 +C2*DP112+C3*DP1 1 3+C4  *DP 114 +C 5* DP  115 +C 6* DP  11 


*  1  6  +  C7  *DP  1 1  7  <■ 
BP-BP  -  AR  *  (c’y. 
*G(2,11)  *  CP  3 ! 
*, 8) *SP5-G(4, 
*■7)  *  (GUI,  7)  * 
' P 1 1 9  +  FM ( 3 )  *  ( 
>9,  11)  ‘CPU  ,  * 


C3*DP118‘-C9‘DP119+C10*DP1110+C11*DPU11)  1  • 

(2)  *  (GUI, 2)  *SP2-G(1,  11)  *CP2 )  *PU2+FM(3)  *  (G(U,3)  ‘SP3-1 
* P 1 1 3  +  FM ( 4 ) * (G ( 1 1 , 4 )  * SP4 -G ( 3,  1 1 ) *CP 4 ) * P 1 1 4+FM  ( 5 ) *  (G ( 1 1 1 
1  I)  *  CP  3 )  *P  1 1  5+FM  (  6)  MG  (11, 6)  *SP6-G(S,  11)  *CP6)  *PU6+FM(1 
Si  7 - G <  U  11)  *CP7) *P1 17+FM (8) * (G ( 1 1 , 8 ) * SP 8-G ( 7 , 1 1 ) *CP8) *1- 
idi,  9)  *Gi  3-0(8,  11)  *CF9)  *P  1 1 9  +  FM  (10)  *  (G<11,  10)  ‘SP10-GU 
P  U  !  :■  *  PM  (1  1 )  *  (G  ( 1 1 ,  1 1 )  *SP1 1  -G  ( 1 0,  1 1 )  ‘CPI  1 )  *P  1 1 1 1 )  1  ■ 


■00244 

-00245 

-00246 

*00247 

-00248 

*00249 

-002S0 

-00251 

-00252 

-00253 

-00254 

-00255 

-00256 

-00257 

-00258 


■00260 

•00261 

-00262 

-00263 

■00264 

-00265 

-00266 

-00267 

-00268 

-00269 

-00270 

-00271 

-00272 

-00273 

-00274 

-00275 

-00276 

-00277 

•00278 

-00279 

-00280 

-00281 

-00282 

-00283 

-00284 

-00285 

•00286 

•00287 

•00288 

•00289 

•00290 

-00291 

■00292 

■00293 

■00294 

■00295 

•00296 

•00297 

■00298 

■00299 

■00300 

■00301 

■00302 

•00303 

•00304 

■00305 
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IF  (NMAX.LE.il)  RETURN 
WRITE  (*, 2)  NMAX 
RETURN 
C 

2  FORMAT (57H0  ERROR,  THIS  MODEL  ONLY  FOR  NMAX=<11,  CALL  WAS  FOR  NMAX 
*=  15) 

C  MAKE  A  SMOOTH  FIT  BETWEEN  TRUNCATED  TERMS. 

10  CON= (R-CONA (N) ) / (CONA(N-l) -CONA(N) ) 

11  BR=BR+ (XR-BR) *CON 
BT=BT+ (XT-BT) *CON 

BP= (BP+ (XP-BP) *CON) /ST 
RETURN 

20  CON= (R-CONA (N) ) / (CONA(N-l) -CONA(N) ) 

2 1  BR=XR+ ( BR-XR) *CON 
BT=XT+ (BT-XT) *CON 

BP= (XP+ (BP-XP) *CON) /ST 
RETURN 

END  1-00480 


1-00306 

1-00476 

1-00477 
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SUBROUTINE  FLDCOF (YEAR, DIMO, NMAXI ) 


C 

c 

c 

n 

c 

c 

c 


c 

c 


DETERMINES  COEFFICIENTS  AND  DIPOL  MOMENT  FROM  IGRF  MODELS 

INPUT:  YEAR  DECIMAL  YEAR  FOR  WHICH  GEOMAGNETIC  FIELD  IS  TO 

BE  CALCULATED 

OUTPUT:  DIMO  GEOMAGNETIC  DIPOL  MOMENT  IN  GAUSS  (NORMALIZED 
TO  EARTH'S  RADIUS)  AT  THE  TIME  (YEAR) 

THIS  ROUNTINE  WAS  INITIALLY  WRITTEN  BY 

D.  BILITZA,  NSSDC,  GSFC,  CODE  633,  GREENBELT,  MD  20771, 

(301)286-9536  NOV  1987. 

MODIFIED  BY  K.  A.  PFITZER  MDSSC  TO  WORK  WITH  GAUSS  NORMALIZED  COEFF . 


C 


c 

C-- 


c 


101 


CHARACTER* 11 
DIMENSION 


DOUBLE  PRECISION 
COMMON /MODEL/ 
COMMON/GENER/ 
DATA 


DATA 


DATA 


FILMOD,  FIL1 ,  FIL2 
GH 1(11,11),  GH2 (11,  11)  , 

DTEMOD (10) , FILMOD(IO) 

X, F0 , F 
G ( 1 1 , 11) 

UMR, ERAD , AQUAD , BQUAD 
/ 'dgrf45 .dat ' ,  'dgrf50.dat*, 

'dgrf60.dat',  'dgrf65.dat', 
'dgrf75.dat',  'dgrf80.dat', 
'  igrf 85s .dat ' / 

,  1950.,  1955.,  I960., 
1975.,  1980.,  1985.,  1990./ 


F I LMOD 
'dgrf55.dat 
'dgrf70 .dat 
' igrf 85 .dat 1 , 
DTEMOD  /  1945 
1965.,  1970., 
LOLD/0/ 


IU  =  11 

DETERMINE  IGRF-YEARS  FOR  INPUT-YEAR 
TIME  =  YEAR 
IYEA  =  INT (YEAR/5 . ) *5 
L  =  (IYEA  -  1945) /5  +  1 


IF  (L.NE.LCLD)  THEN 
LOLD=L 

IF(L.LT.l)  L- 1 
IF ( L . GT . 9 )  L=9 
DTE  1  =  DTEMOD (L) 

FIL1  =  FILMOD (L) 

DTE2  =  DTEMOD (L+l) 

FIL2  =  F ILMOD ( L+ 1 ) 

GET  IGRF  COEFFICIENTS  FOR  THE  BOUNDARY  YEARS 
CALL  GETGAU  (IU,  FIL1,  NMAXI,  ERAD,  GH1,  IER) 

IF  (IER  .NF..  0)  THEN 

WRITE  (*,101)  IU, FIL1, NMAXI, ERAD, IER 

FORMAT  (//'  Error  in  subroutine  FELDCOF'/ 

1  '  IU,  F 1 L 1 ,  NMAXI,  ERAD,  IER:'/ 

2  110, All, 110, 1PE12.3, 110) 

STOP 

ENDIF 

CALL  GETGAU  (IU,  FIL2,  NMAX2 ,  ERAD,  GH2,  IER) 

IF  (IER  .NE.  0)  THEN 

WRITE  (*,102)  IU,FIL2,NMAX2, ERAD, IER 

FORMAT  (//'  Error  in  subroutine  FELDCOF'/ 

1  '  IU,  FIL2,  NMAX2 ,  ERAD,  IER:'/ 

2  110, All, 110, 1PE12.3, 110) 

STOP 

ENDIF 

ENDIF 

T ETEPMINE  IGPF  COEFFICIENTS  FOR  YEAR 
IF  (L  .  LE.  9)  THEN- 

CALL  C INTER  (YEAR,  DTE  1 ,  NMAXI,  GH1 ,  DTE2 , 

1  NMAX2,  OH2,  NMAXI,  G) 

ELSE 
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CALL  EXTRAP  (YEAR,  DTE1,  NMAX1,  GH1,  NMAX2, 
1  GH2,  NMAXI,  G) 

ENDIF 

C —  DETERMINE  MAGNETIC  DIPOL  MOMENT 

F0=G (2, 1) **2+G (2,2) **2+G(l,2) **2 
DIMO=SQRT(FO) 

RETURN 

END 
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SUBROUTINE  GETGAU  < IU,  FSPEC,  NMAX,  ERAD,  G,  IER) 


C 

C  Reads  spherical  harmonic  coefficients  from  the  specified 

C  file  into  an  array  and  converts  the  coefficients  to  Gauss 

C  normalized  coefficients. 


C 

c 

r 

r 

C 


r 

C 

n 

c 

c 

c 


Input : 

IU 

FSPEC 

Output : 
NMAX 
ERAD 


GH 


IER 


-  Logical  unit  number 

-  File  specification 


-  Maximum  degree  and  order  of  model 

-  Earth's  radius  associated  with  the  spherical 
harmonic  coefficients,  in  the  same  units  as 
elevat ion 

-  Gauss  quasi-normal  internal  spherical 
harmonic  coefficients 

-  Error  number:  =  0,  no  error 

=  -2,  records  out  of  order 
=  FORTRAN  run-time  error  number 


CHARACTER  FSPEC* (*) 
DIMENSION  G ( 1 1 , 1 1 ) 


Open  coefficient  file.  Read  past  first  header  record. 
Read  degree  and  order  of  model  and  Earth's  radius. 


OPEN  (IU,  FI LE=FSPEC,  STATUS® 'OLD ' ,  IOSTAT=IER,  ERR=999) 
C  1  READONLY) 

DO  10  1=1,11 
DO  10  J= 1 ,  1  1 
10  G ( I, J) =0 . 

READ  (IU,  *,  IOSTAT = IER,  ERR=999) 

READ  (IU,  *,  IOSTAT® IER,  ERR=999)  MAXN,  ERAD 

IF (MAXN . GT . 10) MAXN® 10 
DO  30  NN= 1 , MAXN 
DO  2  0  MM® 0 , NN 

READ  (IU,  *,  IOSTAT  IER,  ERR=999)  LN,  LM,  GNM,  I1NM 
IF  (NN  .  NE  .  I.N  .  MM  .  NE  .  LM)  THEN 
IER=-2 

_?c  ro  9  '  .■* 

ENDIF 

N=LN+ 1 
M®LMi 1 
G  (N,M) =GNM 
IF(LM.EQ.C)  goto  20 
G  (I,M,  N)  -  HNM 

m.ntinuf. 

3 ;  CONTINUE 

NMAX  =  MAXN  ♦  '. 

:.  v>  ■  r  t  to  Ii:-  •  :  :t  1 1  i  r « •  i 

DO  Lb  N -  1 ,  ' •  M A X 
5  5  M--1  N'"AX 

CALL  CONVPT  t  r. ,  -  )  ,  r,  f  r.,  1  ) 
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55 

CONTINUE 

999 

CLOSE  (IU) 

RETURN 

END 

SUBROUTINE  CONVRT (G, I , L, K) 

DIMENSION  Sill,  ID 
LOGICAL  NEXT 
DATA  NEXT/ .FALSE. / 

IF  (NEXT)  GOTO  2 
NEXT= . TRUE . 

S  ( 1 , 1 ) =-l . 

DO  I  N  =  2, 1  I 

S (N, I) =S (N-l , 1 ) 'FLOAT (2*N-3) /FLOAT (N-l) 

S  (1,N) =0. 

J=2 

DO  1  M=2 , N 

S (N, M) =S (N, M-l ) *SQRT ( (FLOAT (N-M+l) *J) /FLOAT (N+M-2) ) 
S (M-l , N) =S ( N , M ) 

J=1 

IF(K.GT.l)  GOTO  3 
G=G*S (I, L) 

RETURN 
G=G/S (I, L) 

RETURN 

END 
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SUBROUTINE  CINTRP  (DATE,  DTE1,  NMAX1,  GH1,  DTE2 , 
1  NMAX2,  GH2,  NMAX,  GH) 


C 

C  Interpolates  linearly,  in  time,  between  two  spherical 
C  harmonic  models. 

C 

C  Input : 

C  DATE 

C  DTE1 

C  NMAX1 

C  GH1 

C 

C  DTE2 

C  NMAX 2 

C  GH2 

C 
C 

C  Output : 

C  GH 

C  NMAX 

C 

C  ============= 

DIMENSION  GH 1(11, 11),  GH2(11,11),  GH(11,11> 

NMAX=MAX0 (NMAX1 , NMAX2 ) 

FACTOR= (DATE-DTE1 ) / (DTE2-DTE1 ) 

DO  234  J=  1, 11 
DO  234  1=1,  11 

234  GH (I, J)  =  GH1 ( I , J)  +  FACTOR  *  <GH2(I,J)  -  GH1(I,J)) 

RETURN 
END 


-  Date  of  resulting  model  (in  decimal  year) 

-  Date  of  earlier  model 

-  Maximum  degree  and  order  of  earlier  model 

-  Gauss  quasi-normal  internal  spherical 
harmonic  coefficients  of  earlier  model 

-  Date  of  later  model 

-  Maximum  degree  and  order  of  later  model 

-  Gauss  quasi-normal  internal  spherical 
harmonic  coefficients  of  later  model 


-  Coefficients  of  resulting  model 

-  Maximum  degree  and  order  of  resulting  model 
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SUBROUTINE  EXTRAP  (DATE,  DTE1 ,  NMAX1,  GH1 ,  NMAX2, 
1  GH2 ,  NMAX,  GH) 


C 

C 

c 

c 

c 


r 

n 

r 

C 

r 

C 

c 

c 


Extrapolates  linearly  a  spherical  harmonic  model  with  a 
rate-of-change  model. 


Input : 

DATE 
DTE  1 
NMAX1 
GH1 

NMAX2 

GH2 


Date  of  resulting  model  (in  decimal  year) 

Date  of  base  model 

Maximum  degree  and  order  of  base  model 
Gauss  quasi-normal  internal  spherical 
harmonic  coefficients  of  base  model 
Maximum  degree  and  order  of  rate-of-change 
model 

Gauss  quasi-normal  internal  spherical 
harmonic  coefficients  of  rate-of-change  model 


Output : 

GH  -  Coefficients  of  resulting  model 

NMAX  -  Maximum  degree  and  order  of  resulting  model 


DIMENSION  GH1 (11, 11) ,  GH2(11,11),  GH(11,11) 

NMAX=MAX0 ( NMAX 1 , NMAX2 ) 

FACTOR  =  (DATE  -  DTE1 ) 

DO  567  J=l, 11 
DO  567  I  =  1,11 

567  GH ( I, J)  =  GH 1 ( I , J)  +  FACTOR  *  GH2(I,J) 

RETURN 

END 
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Appendix  B 

External  Magnetic  Field  Routines 

These  listings  are  the  1 977  Olson-Pfitzer  tilt  dependent  routine  including  the  routine  that 
combines  the  internal  and  the  external  field.  These  routines  are  included  here  for 
completeness.  They  were  not  developed  or  changed  as  a  part  of  this  effort.  They  are, 
however,  necessary  in  order  for  the  remaining  software  to  function  properly. 
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n  o  o  n  o  o  o  an  n  no  oonoononnnooo  non 


BMNIGRF  --  A  TEST  ROUTINE  TO  CHECK  THE  OPERATION  OF  THE  TILT 
DEPENDENT  MODEL  AND  IT  COMBINATION  WITH  THE  IGRF  MAIN  FIELD 

DIMENSION  X (3) , B (3) 

COMMON /BXYZCM/ YEAR, DAYYR, UT, KODE, JSW 


SET  UP  THE  YEAR  FOR  THE  MAIN  FIELD  ROUTINE 
YEAR=1 98  5 . 

SET  THE  SWITCH  TO  USE  EXTERNAL  PLUS  INTERNAL  FIELD 


JSW=1 


SET  THE  SWITCH  TO  USE  INPUT  AND  OUTPUT  IN  CARTESIAN  COORDS 
KODE= 1 

SET  UP  A  CARTESIAN  COORDINATE  TEST  LOOP 

SET  UP  DATE  AND  TIME 
DO  200  ID=1,2 
DAYYR=90*ID 
DO  190  IUT=1 ,  3 
UT=IUT  *6-6 


PRINT  PAGE  HEADER 
WRITE  (6,110) 

11C  FORMAT (77H1DAYOFYR  UT  X  Y  Z  BX  BY 

*  BZ  BMAG, / ) 

SET  UP  POSITION  IN  CARTESIAN  COORDS 
DO  180  IZ=1, 3 
X ( 3 ) =3  * IZ- 6 
DO  170  IY=1 , 3 
X (2) =3*IY-6 
DO  160  IX=1, 6 
X  ( 1 ) =4 * IX- 1 4 

GET  THE  MAGNETIC  FIELD  VALUES 
CALL  BMNEXT (X, B, BMAG) 

WRITE  (6,120)  DAYYR, UT,X,B, BMAG 
120  FORMAT (F6 . 0,  4F8 .2, 4F10. 5) 

160  CONTINUE 
170  CONTINUE 
180  CONTINUE 
190  CONTINUE 
200  CONTINUE 

SET  UP  FOR  SPHERICAL  COORDINATES 

000  KODE=2 

SET  DATE  AND  TIME 
DO  300  ID= 1 , 2 
DAYYR=90  *  ID 
DO  290  IUT- 1 , 3 
UT=IUT*6-6 


PRINT  PAGE  HEADER 
WRITE  (6,210) 


FORMAT (77H1DAYOFYR  UT 

R 

THETA 

PHI 

BR 

BTHE 

TA  BPHI  BMAG,/) 
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C  SET  UP  POSITIONS  IN  SPHERICAL  COORDS 

DO  280  IR-1, 3 
X ( 1) =IR*3 
DO  270  IT=1 , 3 
X  (2) =IT*45 
DO  260  IP-1,  6 
X (3) - (IP-1) *60 

GET  THE  MAGNETIC  FIELD 
CALL  BMNEXT (X, B, BMAG) 

WRITE (6, 120)  DAYYR, UT, X, B, BMAG 
260  CONTINUE 
270  CONTINUE 
280  CONTINUE 
290  CONTINUE 
300  CONTINUE 
END 
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SUBROUTINE  BMNEXT (XX, B, BMAG) 

C 

C  PURPOSE 

C  TO  DETERMINE  THE  MAIN  MAGNETIC  FIELD  PLUS  THE  EXTERNAL 

C  FIELD 

C 

C  METHOD 

C  DETERMINES  THE  VECTOR  MAGNETIC  FIELD  IN  GEOGRAPHIC 

C  COORDINATES  USING  A  SPHERICAL  COORDINATE  EXPANSION  OF  THE 

C  EARTHS  INTERNAL  FIELD  AND  A  CARTESIAN  COORDINATE  EXPANSION 

C  OF  THE  BOUNDARY,  TAIL  AND  RING  CURRENT  FIELDS  IN  SOLAR 

C  MAGNETIC  COORDINATES 


C 

r 

C 

c 

c 


r 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r* 

c 

c 

c 

c 

c 

c 

c 

r 

r 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


INPUT  --  ARGUMENT  LIST 

XX  A  REAL  ARRAY  CONTAINING  THE  POSITION  IN  GEOGRAPHIC 

COORDINATES 
IF  KODE  =  1 

XX(1)=X,  XX (2 ) =  Y,  XX (3) =  Z,  WHERE  X,  Y,  Z  ARE  IN  EARTH 
RADII.  THE  DIRECTION  OF  Z  IS  ALONG  THE  EARTHS  ROTATION 
AXIS  TOWARDS  THE  GEOGRAPHIC  NORTH  POLE.  THE  DIRECTION 
OF  X  IS  TO  THE  GREENWHICH  MERIDIAN  IN  THE  EQUATORIAL 
PLANE.  THE  Y  AXIS  IS  IN  THE  EQUATORIAL  PLANE  NORMAL 
TO  X  AND  Z  IN  A  RIGHT  HANDED  SENSE. 

IF  KODE  =  2 

XX ( I ) =R,  GEOCENTRIC  RADIUS  IN  EARTH  RADII, 

XX (2) =THETAG,  COLATITUDE  IN  DEGREES, 

XX (3) =PHIG,  LONGITUDE  IN  DEGREES 

INPUT  --  COMMON  BLOCK  BXYZCM 

UT  THE  CURRENT  UNIVERSAL  TIME  IN  HOURS 

KODE  A  FLOW  CONTROL  VARIABLE.  KODE  EQUAL  TO  ONE  MEANS  THAT 
INPUT  AND  OUTPUT  ARE  IN  CARTESIAN  COORDINATES.  KODE 
EQUAL  TO  TWO  MEANS  THAT  INPUT  AND  OUTPUT  ARE  SPHERICAL 
COORDINATES . 

DAYYR  THE  NUMBER  OF  THE  DAY  OF  YEAR 

JSW  A  FLOW  CONTROL  VARIABLE.  IF  JSW  IS  LESS  THAN  ZERO,  THE 
THE  FIELD  IS  COMPUTED  USING  THE  INTERNAL  FIELD  ONLY. 

IF  JSW  IS  GREATER  THAN  OR  EQUAL  TO  ZERO  THE  FIELD 
WILL  BE  COMPUTED  USING  THE  INTERNAL  PLUS  EXTERNAL 
FIELD . 

YEAR  THE  YEAR  USED  BY  THE  INTERNAL  MAGNETIC  FIELD  ROUTINE 
TO  TAKE  INTO  ACCOUNT  THE  SECULAR  VARIATIONS 
(E.G.  JULY  15,  1964  =  1964.54) 

NOTE****  YEAR  SHOULD  BE  CHANGED  ONLY  EVERY  FEW  DAYS  OR 
MONTHS.  NEW  FIELD  COEFFICIENTS  MUST  BE  COMPUTED  FOR 
EVERY  CHANGE  IN  YEAR.  THIS  COULD  CAUSE  A  LARGE  INCREASE 
IN  COMPUTER  TIME.  THE  EARTHS  FIELD  CHANGES  ONLY  ABOUT 
.001  GAUSS/YEAR  AT  THE  EARTHS  SURFACE. 

OUTPUT  --  ARGUMENT  LIST 

B  A  REAL  ARRAY  CONTAINING  THE  COMPONENTS  OF  THE  MAGNETIC 

FIELD  IN  GAUSS  AT  THE  CURRENT  POSITION  AND  TIME 
IF  KODE  =  1 

Li  ( 1 )  =EX,  B  (2 )  =BY,  B  (3)  =BZ  THE  CARTESIAN  COMPONENTS 
OF  THE  MAGNETIC  FIELD  IN  GEOGRAPHIC  COORDINATES 
IF  KODE  =  2 

B ( 1 ) =BR,  RADIAL  COMPONENT  OF  THE  FIELD,  POSITIVE  IN  THE 
DIRECTION  OF  INCREASING  RADIUS. 

B (2 ) -BTHETA,  COMPONENT  IN  LATITUDE,  POSITIVE  IN  THE 
DIRECTION  OF  INCREASING  COLATITUDE 

B  (  3 ) =BPH  I ,  COMPONENT  IN  LONGITUDE,  POSITIVE  IN  THE 
DIRECTION  OF  INCREASING  LONGITUDE. 

BMAG  THE  MAGNITUDE  OF  THE  MAGNETIC  FIELD  VECTOR  IN  UNITS  OF 

GA'-SS  . 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


OUTPUT  —  COMMON  BLOCK  BXYZCM 

XMLAT  THE  MAGNETIC  LATITUDE  AT  THE  CURRENT  POSITION  IN  RADIANS 

SUBROUTINE  CONSTANTS 

PICON  THE  NUMBER  OF  DEGREES  PER  RADIAN 

SIN  D  THE  SINE  OF  THE  COLATITUDE  OF  THE  DIPOLE  AXIS 

COS  D  THE  COSINE  OF  THE  COLATITUDE  OF  THE  DIPOLE  AXIS 

C69  COSINE  OF  69 

S69  SINE  OF  69 

CALLING  SUBROUTINES 
SUBROUTINE  INVARM 

SUBROUTINES  REQUIRED 
SUBROUTINE  BXYZMU 
SUBROUTINE  ANGLE 
SUBROUTINE  SPIGRF 

VARIABLES 

AOR  INVERSE  OF  RADIUS  VECTOR  (AOR=l . /R) 

BGMX , BGMY , BGMZ  INTERMEDIATE  VALUES  OF  THE  MAGNETIC  FIELD 
VECTOR  DURING  COORDINATE  TRANSFORMATION 
BMX, BMY, BMZ  EXTERNAL  MAGNETIC  FIELD  IN  GEOMAGNETIC  COORDINATES 
BP, BR, BT  COMPONENTS  OF  INTERNAL  FIELD  IN  SPHERICAL  COORDINATES 
BP  IS  LONGITUDINAL  COMPONENT 
BR  IS  RADIAL  COMPONENT 
BT  IS  LATITUDINAL  COMPONENT 

BX, BY, BZ  CARTESIAN  COMPONENT  OF  EXTERNAL  FIELD  IN  GEOGRAPHIC 
COORDINATES 

CP  COSINE  OF  COLATITUDE 

CPS  COSINE  OF  HOUR  ANGLE  TO  GET  FROM  SOLAR  MAGNETIC  TO 
GEOMAGNETIC  COORDINATES 
CT  COSINE  OF  GEOGRAPHIC  LONGITUDE 

DAYLST  LAST  DAY  FOR  WHICH  TILT  AND  HOUR  ANGLE  WERE  UPDATED 
NMAX  MAXIMUM  NUMBER  OF  TERMS  USLb  bY  INTERNAL  FIELD  ROUTINE 
SET  UP  BY  INTERNAL  FIELD  ROUTINE 
PHIG  GEOGRAPHIC  COLATITUDE 
R  RADIUS  VECTOR  TO  POSITION  POINT 

R2  R*  *2 

SP  SINE  OF  COLATITUDE 

SPS  SINE  OF  HOUR  ANGLE  TO  GET  FROM  SOLAR  MAGNETIC  TO 

GEOMAGNETIC  COORDINATES 
ST  SINE  OF  LONGITUDE 

THETAG  GEOGRAPHIC  LONGITUDE 
TILT  TILT  OF  THE  DIPOLE  AXIS 

UTLST  LAST  UNIVERSAL  TIME  FOR  WHICH  TILT  AND  HOUR  ANGLE  WERE 
UPDATED 

X  A  REAL  ARRAY  HOLDING  THE  POSITION  VECTOR  IN  SOLAR 

MAGNETIC  COORDINATES 

XP, YP, ZP  POSITION  VECTOR  IN  GEOMAGNETIC  COORDINTES 
XPP, YPP  INTERMEDIATE  POSITION  COMPONENT  DURING  COORDINATE 
TRANSFORMATION 

YEARI  TRANSMITS  THE  YEAR  TO  THE  INTERNAL  FIELD  ROUTINE 
VERSION  10/25/77 

FOR  MORE  INFORMATION  CALL  OR  WRITE  K.  A.  PFITZER  AT  MCDONNELL 
DOUGLAS  ASTRONAUTICS  CO.  5301  BOLSA  AVE,  HUNTINGTON  BEACH  CALIF. 
PHONE  (714)  896-3231. 

DIMENSION  X(3) ,B(3) ,XX(3) 

COMMON/ BXYZCM/ YEAR, DAYYR, UT, KODE, JSW 

COMMON  /GCOM/  ST, CT, SP, CP, AOR, BT, BP, BR, NMAX, YEARI 

DATA  PICON/57.29577951/, SIND, COSD/. 2027872954, . 9792228106/, 
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*S69,C69/. 933580 42 65,  . 35836794 95/ , UTLST, DAYLST/2 *  12 3456  . / 

UPDATE  THE  ROTATION  HOUR  ANGLE  AND  TILT  ANGLE  IF  THE  UNIVERSAL 
TIME  OR  DAY  OF  YEAR  HAS  CHANGED  SINCE  THE  LAST  CALL 

IF (UT.EQ. UTLST. AND. DAYYR . EQ . DAYLST )  GO  TO  1 

UTLST=UT 

DAYLST=DAYYR 

CALL  ANGLE  (TILT, SPS, CPS) 

1  IF(KODE.GT.l)  GO  TO  3 

DETERMINE  THE  SPHERICAL  COORDINATES  OF  POSITION  IF  CARTESIAN 
COORDINATES  WERE  ENTERED 

X(1)-XX(1) 

X (2) =XX  (2) 

X (3) =XX (3) 

R2  =X ( 1 ) *  * 2  +  X  ( 2 ) *  *  2 
R=SQRT ( X ( 3 ) **2+R2) 

R2=SQRT (R2 ) 

CT=X (3) /R 
ST-R2/R 
CP=X<1) /R2 
SP=X  (2) /R2 
GO  TO  5 

DETERMINE  THE  CARTESIAN  COORDINATES  OF  POSITION  IF  SPHERICAL 
COORDINATES  WERE  ENTERED 

3  R=XX  ( 1 ) 

THETAG=XX (2) /PICON 
PHIG=XX (3) /PICON 
CT=COS (THETAG) 

ST=SIN (THETAG) 

CP=COS (PHIG) 

SP=SIN (PHIG) 

X(l) =R*ST*CP 
X (2) =R*ST*SP 
X ( 3 ) =R*CT 
5  BX=0 . 

BY=0  . 

BZ  =  0  . 

IF  THE  EXTERNAL  MAGNETIC  FIELD  IS  TO  BE  USED  IN  THE  COMPUTATION, 
COMPUTE  THE  SOLAR  MAGNETIC  COORDINATES 

IF(JSW.LT.O)  GO  TO  9 

FIRST  ROTATION  IS  ABOUT  THE  Z-AXIS  TROUGH  AN  ANGLE  OF  291  DEGREES 
(THE  LONGITUDE  OF  THE  MAGNETIC  NORTH  POLE) 

XPP=X(1) *C  6  9-X  ( 2 ) *S69 
YPP=  X ( 1 >  'S69+X  (2)  *C69 

SECOND  POTATION  IS  ABOUT  THE  NEW  Y-AXIS  THROUGH  AN  ANGLE  OF  11.7 
DEGREES  (THE  COLATITUDE  OF  THE  MAGNETIC  NORTH  POLE) 

ZP  =  XPP  *S IND  +  X ( 3 )  ‘COSD 
XP=XPP ‘COSD-X ( 3 ) *  S IND 
YP= YPP 

C  ROTATION  IS  ABOUT  THE  MAGNETIC  Z-AXIS  THROUGH  THE  HOUR  ANGLE  OF 

C  THE  SUN  FROM  THF.  PRIME  MAGNETIC  MERIDIAN  (NEGATIVE  ROTATION) 
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X ( 1 ) -XP*CPS-YP*SPS 
X(2)«XP*SPS+YP*CPS 
X ( 3) =ZP 

DETERMINE  THE  EXTERNAL  MAGNETIC  FIELD  USING  A  TILT  DEPENDENT 
MAGNETIC  FIELD 

CALL  BXYZMU (X, B, TILT) 

THE  CARTESIAN  COMPONENTS  OF  THE  FIELD  ARE  IN  SOLAR  MAGNETIC 
COORDINATES.  THE  COMPONENTS  ARE  NEEDED  IN  THE  GEOGRAPHIC 
COORDINATE  SYSTEM 

FIRST  ROTATION  IS  ABOUT  THE  MAGNETIC  Z-AXIS  THROUGH  THE  HOUR 
ANGLE  OF  THE  SUN  TO  THE  PRIME  MAGNETIC  MERIDIAN 
(POSITIVE  ROTATION)  PUTS  RESULTS  INTO  GEOMAGNETIC  COORDINATES 

BMX=B ( 1 ) *CPS+B (2) *SPS 
BMY=~B ( 1 ) *SPS+B (2 ) *CPS 
BMZ=B (3) 

SECOND  ROTATION  IS  ABOUT  THE  MAGNETIC  Y-AXIS  THOUGH  -11.7  DEGREES 
COLATITUDE 

BGMX=BMX*COSD+BMZ*SIND 

BGMY=BMY 

BGMZ=-BMX*SIND+BMZ*COSD 

THIRD  ROTATION  IS  ABOUT  THE  NEW  Z-AXIS  THROUGH  -291  DEGREES 

BX=BGMX*C69+BGMY*S69 

BY=-BGMX*S69+BGMY*C69 

BZ=BGMZ 

DETERMINE  THE  MAIN  FIELD 

9  CONTINUE 
AOR=l . /R 
YEARI=YEAR 
CALL  SPIGRF 
IF(KODE.GT.l)  GO  TO  10 

IF  THE  OUTPUT  IS  TO  BE  IN  CARTESIAN  GEOGRAPHIC  COORDINATES  CONVERT 
THE  MAIN  MAGNETIC  FIELD  AND  ADD 

B ( 1 ) = (BX+CP* (ST*BR+CT*BT) -SP*BP) *0.00001 
B (2) = (BY+SP* (ST*BR+CT*BT) +CP*BP) *0.00001 
B (3) = (BZ+CT*BR~ST*BT) *0.00001 
GO  TO  20 

IF  OUTPUT  IS  TO  BE  IN  SPHERICAL  GEOGRAPHIC  CONVERT  THE  EXTERNAL 
FIELD  AND  ADD 

10  B (1) = (BR+ (BX*CP+BY*SP) *ST+BZ*CT) *0.00001 
B (2) = (BT+ (BX*CP+BY*SP) *CT-BZ*ST) *0.00001 
B(3)=(BP+BY*CP-BX*SP) *0.00001 

DETERMINE  THE  MAGNITUDE  OF  THE  FIELD  VECTOR 

20  BMAG=SQRT (B(l)**2+B(2)**2+B(3)**2) 

RETURN 

END 
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SUBROUTINE  ANGLE (TILT, SINPHE, COSPHE) 

C 

C  PURPOSE 

C  THIS  ROUTINE  CALCULATES  THE  ANGLE  BETWEEN  THE  MAGNETIC  DIPOLE 

C  AXIS  AND  THE  SUN-EARTH  LINE  AS  WELL  AS  THE  ROTATION  SINES 

C  AND  COSINES  TO  CONVERT  FROM  GEOMAGNETIC  TO  SOLAR  MAGNETIC 

C  COORDINATES 

C  METHOD 

C  MAGNETIC  COORDINATES  HAVE  THEIR  ORIGIN  AT  THE  CENTER  OF  THE 

C  EARTH  WITH  THE  Z  AXIS  ALLIGNED  THROUGH  THE  GEOMAGNETIC  NORTH 

C  POLE.  IN  GEOMAGNETIC  COORDINATES  THE  X  AXIS  IS  IN  THE 

C  PLANE  PASSING  THROUGH  THE  DIPOLE  AXIS  AND  THE  GEOGRAPHIC 

C  AXIS  (ABOUT  69  DEGREES  WEST  LONG.).  IN  SOLAR  MAGNETIC 

C  COORDINATES  X  AXIS  LIES  IN  THE  PLANE  CONTAINING  THE  SUN 

C  EARTH  LINE  AND  THE  Z  AXIS  (POSITIVE  X  AXIS  HAS  A  LARGE 

C  COMPONENT  IN  THE  SOLAR  DIRECTION) .  THE  Y  AXIS  IS  ORTHOGONAL 

C  TO  THE  X  AND  Z  AXIS  SUCH  THAT  X,  Y  AND  Z  FORM  A  RIGHT 

C  HANDED  SYSTEM.  THE  ECCLIPTIC  COORDINATE  SYSTEM  HAS  ITS 

C  Z  AXIS  ALONG  THE  ECCLIPTIC  NORTH  POLE  (THROUGH  THE  CENTER 

C  OF  THE  EARTH  AND  PERPENDICULAR  TO  THE  EARTHS  ORBITAL  PLANE) 

C  THE  X  AXIS  POINTS  TOWARD  THE  SUN  AND  Y  FORMS  A  RIGHT  HANDED 

C  COORDINATE  SYSTEM.  IN  THIS  ROUTINE  IN  ORDER  TO  REDUCE 

C  COMPUTER  TIME  THE  APPROXIMATION  OF  A  CIRCULAR  EARTH  ORBIT 

C  AROUND  THE  SUN  IS  MADE. 

C 

C  INPUT  --  COMMON  BLOCK  BXYZCM 

C  DAYYR  IS  THE  DAY  OF  YEAR  (1.-366.).  IT  MUST  BE  A  WHOLE 

C  NUMBER.  DAY  1  IS  JANUARY  1. 

C  UT  THE  UNIVERSAL  TIME  IN  HOURS  (0.0000-2-1.00000) 

C  OUTPUT  --  PARAMETER  LIST 

TILT  THE  TILT  OF  THE  DIPOLE  AXIS  IN  DEGREES. 

TILT  =  90.  -  PSI,  WHERE  PSI  IS  THE  ANGLE  BETWEEN 
THE  MAGNETIC  DIPOLE  AXIS  AND  THE  SOLAR  DIRECTION. 
SINPHE  THE  SINE  OF  THE  ROTATION  ANGLE  ABOUT  THE  MAGNETIC 

Z  AXIS  TO  CONVERT  FROM  GEOMAGNETIC  TO  SOLAR  MAGNETIC 
COORDINATES 

COSPHE  THE  COSINE  OF  THE  ROTATION  ANGLE  ABOUT  THE  MAGNETIC 
Z  AXIS  TO  CONVERT  FROM  GEOMAGNETIC  TO  SOLAR  MAGNETIC 
C  COORDINATES 


c 

CONSTANTS 

c 

P 12 

PI  /  2 

f* 

CON 

180.  / 

PI  CONVERTS  RADIANS  TO  DEGREES 

c 

SALF 

SINE  ( 

11.1)  INCLINATION  OF  MAGNETIC  Z  TO  GEOGRAPHIC 

c 

CALF 

COSINE 

(11.7) 

c 

SGAM 

SIN  (2 

3.5)  INCLINATION  OF  ROTATION  AXIS  TO  ECLIPTIC  Z 

c 

CGAM 

COSINE 

(23.5) 

c 

SALF  ' 

3GAH 

c 

SACG 

SALF  * 

CGAM 

r 

CASG 

CALF  * 

SGAM 

c 

CACG 

CALF  * 

CGAM 

r* 

W 

EARTHS 

ANGUI-AR  ROTATION  FREQUENCY  CORRECTED  FOR  ITS 

c 

ONCE  A 

YEAR  ROTATION  ABOUT  THE  SUN  (UNITS  ARE  1 /HOURS 

C  CALLING  SUBROUTINES 

C  SUBROUTINE  BMNF.XT 

C  VARIABLES 

C  WT  INSTANTANEOUS  ROTATION  ANGLE  AT  THE  SPECIFIED  UNIVERSAL 

C  TIME  AND  DAY  OF  YEAR 

C  CWT  WT/ 366. 256 

C  XX,  XY,  XZ  "  N'ENTS  OF  THE  GEOMAGNEITIC  X  AXIS  IN  ECCLIPTIC 
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C  COORDINATES 

C  ZX, ZY,ZZ  COMPONENTS  OF  THE  DIPOLE  AXIS  IN  ECCLIPTIC  COORDINATES 

C  OSP, SSMLT, CSMLT, SBWT, CBWT, SMLSWT , SMLCWT, CMLSWT, CMLCWT  ARE 

C  SINES  AND  COSINES  AND  THEIR  PRODUCTS  AND  ARE  SET  UP 

C  TO  MINIMIZE  COMPUTER  TIME 

C 

COMMON/BXYZCM/YEAR, DAYYR, UT, KODE, JSW 
DATA  IFIRST/0/ 

C 

C  THE  FIRST  TIME  TROUGH  THE  SUBROUTINE  SET  UP  THE  FIXED  CONSTANTS 

IF ( IFIRST . NE . 0 )  GO  TO  10 
IFIRST=1 

PI2=ATAN2 (0. , -1 . ) /2 . 

CON=90./PI2 
SALF=SIN (11 . 7/CON) 

CALF=COS (11.7 /CON) 

SGAM=SIN (23 . 5/CON) 

CGAMCOS  (23 .5 /CON) 

SASG=SALF*SGAM 
SACG=SALF*CGAM 
CASG=CALF  *  SGAM 
CACG=CALF  *CGAM 
W=PI2/6 . * ( 1 . +1 . /36S.2S6) 

C 

C  MAIN  ENTRY  POINT.  SET  UP  THE  THE  SINES  AND  COSINES  REQUIRED 

C  BY  THE  TRANSFORMATIONS. 

10  WT= (UT-16 . 6+ (DAYYR-172 . ) *24.) *W 
CWT=-WT/365 .256 
SSMLT=SIN (WT) 

CSMLT=COS (WT) 

SBWT=SIN (CWT) 

CBWT=COS (CWT) 

SMLSWT=SSMLT*SBWT 
SMLCWT=SSMLT*CBWT 
CMLS WT=CSMLT  *  SBWT 
CMLCWT-CSMLT*CBWT 
C 

C  DETERMINE  THE  COMPONENTS  OF  THE  DIPOLE  AXIS  IN  ECCLIPTIC 

C  COORDINATES 

ZX=CASG*CBWT+SACG*CMLCWT-SALF* SMLSWT 
Z Y=CASG * SBWT +SACG*CMLSWT+SALF* SMLCWT 
ZZ=CACG-SASG* CSMLT 
C 

C  CALCULATE  THE  TILT  ANGLE 

PSI=ACOS (ZX) 

OSP=l . / (SIN (PSI) ) 

TILT=CON* (PI2-PSI) 

C 

C  DETEMINE  THE  COMPONENTS  OF  THE  GEOMAGNETIC  X  AXIS  IN  ECCLIPTIC 

C  COORDINATES 

XX=C ACG*CMLCWT-SASG* CBWT -CALF  * SMLSWT 
XY=CACG*CMLSWT-SASG  *  SBWT +CALF* SMLCWT 
XZ=-CASG*CSMLT-SACG 
C 

C  OBTAIN  THE  ROTATION  SINES  AND  COSINES 

SINPHE= (XY*  ZZ-XZ  *  ZY) *OSP 

COSPHE= (XX* (ZZ*ZZ+ZY*ZY) -ZX* (XY*ZY+XZ * ZZ ) ) *OSP 

RETURN 

END 
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SUBROUTINE  BXYZMU (XX, BF, TILT) 

VERSION  11/01/76 
PURPOSE 

TO  CALCULATE  THE  CONTRIBUTION  TO  THE  EARTHS  MAGNETIC  FIELD  BY 
SOURCES  EXTERNAL  TO  THE  EARTH.  NO  INTERNAL  FIELD  IS  INCLUDED 
IN  THIS  ROUTINE. 

METHOD 

THE  ROUTINE  INCLUDES  THE  FIELD  CONTRIBUTIONS  FROM  THE 
MAGNETOPAUSE  CURRENTS,  AND  CURRENTS  DISTRIBUTED  THROUGHOUT 
THE  MAGNETOSPHERE  (THE  TAIL  AND  RING  CURRENTS) .  IT  IS  VALID 
FOR  ALL  TILTS  OF  THE  EARTHS  DIPOLE  AXIS  AND  IS  VALID  DURING 
QUIET  MAGNETIC  CONDITIONS. 

A  GENERALIZED  ORTHONORMAL  LEAST  SQUARES  PROGRAM  WAS  USED 
TO  FIT  THE  COEFFICIENTS  OF  A  POWER  SERIES  (INCLUDING 
EXPONENTIAL  TERMS)  THROUGH  FOURTH  ORDER  IN  SPACE  AND 
THIRD  ORDER  IN  TILT.  THIS  EXPANSION  HAS  BEEN  OPTIMIZED 
FOR  THE  NEAR  EARTH  REGION  AND  IS  VALID  TO  15  EARTH  RADII. 
OUTSIDE  OF  THIS  REGION  THE  FIELD  DIVERGES  RAPIDLY  AND  A 
TEMPLATE  SETS  THE  FIELD  TO  ZERO.  IN  ORDER  TO  IMPROVE 
COMPUTATIONAL  SPEED  THE  FIELD  IS  SET  TO  ZERO  BELOW  2  EARTH 
RADII.  (IN  THIS  REGION  THE  EARTHS  INTERNAL  FIELD  DOMINATES 
AND  THE  VARIATIONS  EXCPRESSED  BY  THIS  EXPANSION  IS  NOT 
SUFFICIENTLY  ACCURATE  THE  PREDICT  VARIATIONS  ON  THE  EARTHS 
SURFACE) 

THE  POWER  SERIES  REPRESENTING  THE  MAGNETIC  FIELD  IS 
BX=SUM  OVER  I,J,K  OF  (  A ( I , J, K) *X* * ( 1-1 ) * Y* * (2 * J-2 ) *Z* * (K-l ) 

+  B (I,  J,K) *X** (I-l) * Y*  * (2  * J-2 ) *Z**(K-1) *EXP (-.06*R**2) ) 
I  GOES  FROM  1  TO  5,J  FROM  1  TO  3,  K  FROM  1  TO  5 

THE  SUM  OF  I  +  2*J  +  K  IS  LESS  THAN  OR  EQUAL  TO  9 

BY=SUM  OVER  I,J,K  OF  (  C ( I , J, K) *X* * ( I-l ) * Y* * ( 2  * J-l ) *Z *  * (K-I ) 

+  D ( I , J, K) *X*  * ( I -1 ) *Y** (2* J-l) *Z** (K-l) *EXP (- . 06*R**2) ) 
I  GOES  FROM  1  TO  5,  J  FROM  1  TO  3,  K  FROM  1  TO  5 
THE  SUM  OF  I  +  2*J+1  +  K  IS  LESS  THAN  OR  EQUAL  TO  9 
BZ  =  SUM  OVER  I , J, K  OF  (  E ( I , J, K) *X* * ( I -1 ) * Y* * (2 * J-2 ) * Z* * (K-l ) 

+  F ( I ,  J,K)  *X** (I-l)  *Y*  * (2*J-2 ) *Z** (K-l) *EXP(-.0  6*R**2) ) 
I  GOES  FROM  1  TO  5,  J  FROM  1  TO  3,  K  FROM  1  TO  5 
THE  SUM  OF  I  +  2*J  +  K  IS  LESS  THAN  OR  EQUAL  TO  9 

THE  COEFFICIENTS  A-F  ARE  DEPENDENT  ONLY  ON  POSITION  AND 

ARE  RECALCULATED  EACH  TIME  THE  TILT  OF  THE  DIPOLE  IS  CHANGED. 
THE  COEEF IC IENTS  A-F  ARE  DETERMINED  FROM  THE  TILT  DEPENDENT 
CONSTANTS  AA-FF  BY  THE  FOLLOWING  EXPRESSIONS 
A(I,  J,K)  =  AA (I,  J,  K, 1) *TILT** (K-l- (K-l) /2*2) 

*A/.  ( I,  J,K,  2)  *TILT**  (K+l-  (K-l)  /2*2) 

B ( I, J, K) =BB . 

C(I,  J,K)  =CC . 

D ( I , J, K) =DD . 

E  (I,  J,  K)  =EE  (I,  J,  K,  1)  *TILT**  (K-  (K)  /  2  *  2 ) 

+EE ( I , J, K, 2 ) *TILT*  * (K+2-(K) / 2*2) 

F (I,  J, K) =FF . 


INPUT  --  CALLING  SEQUENCE 

XX  A  REAL  ARRAY  GIVING  THE  POSITION  WHERE  THE  MAGNETIC 

FIELD  IS  TO  BE  EVALUATED.  XX(1),  XX(2),  XX(3)  ARE 
RESPECTIVELY  THE  X,  Y,  AND  Z  SOLAR  MAGNETIC 
COORDINATES  IN  EARTH  RADII.  Z  IS  ALONG  THE  EARTHS 
NORTH  DIPOLE  AXIS,  X  IS  PERPENDICULAR  TO  Z  AND  IN  THE 
PLANE  CONTAINING  THE  Z  AXIS  AND  THE  SUN-EARTH  LINE, 

Y  IS  PERPENDICULAR  TO  X  AND  Z  FORMING  A  RIGHT  HANDED 
COORDINATE  SYSTEM.  X  IS  POSITIVE  IN  THE  SOLAR  DIRECTION. 

TILT  I"  THE  TILT  OF  THE  DIPOLE  AXIS  IN  DEGREES.  IT  IS 
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THE  COMPLEMENT  OF  THE  ANGLE  BETWEEEN  THE  NORTH  DIPOLE 
AXIS  AND  THE  SOLAR  DIRECTION  (PSI).  TILT=90-PSI . 

OUTPUT  —  CALLING  SEQUENCE 

BF  A  REAL  ARRAY  CONTAING  THE  X,  Y,  AND  Z  COMPONENTS  OF 

THE  MAGNETOSPHERIC  MAGNETIC  FIELD  IN  GAMMA.  BF(1), 

BF (2)  AND  BF (3)  ARE  THE  BX,  BY,  BZ  COMPONENTS. 


CONSTANTS 

AA, BB, CC, DD, EE, FF  ARE  REAL  ARRAYS  CONTAING  THE  TILT  DEPENDED 
COEFFICIENTS. 

AA ( I, J, K, L)  ARE  STORED  SUCH  THAT  L  VARIES  MOST  RAPIDLY 
FOLLOWED  IN  ORDER  BY  K,  J  AND  I.  I  VARIES  THE  SLOWEST 
THE  ARRAY  IS  CLOSE  PACKED  AND  ALL  COEFFICIENTS  THAT 
ARE  ZERO  BECAUSE  OF  SYMMETRY  OR  BECAUSE  THE  CROSS  TERM 
POWER  IS  TOO  LARGE  ARE  DELETED. 


VARIABLES 

A, B, C, D, E, F  THE  TILT  INDEPENDENT  COEFFICIENTS.  THEIR  USE 
IS  DESCRIBED  UNDER  METHOD. 

ITA  A  REAL  ARRAY  WHICH  CONTAINS  THE  SYMMETRY  OF  THE  TILT 
DEPENDENCE  FOR  EACH  OF  THE  A  AND  B  COEFFICIENTS 
ITA ( 1 )  HAS  THE  SYMMETRY  INFORMATION  FOR  A(l, 1,1,1) 

AND  A ( 1 , 1, 1,2) 

ITA (2 )  HAS  THE  SYMMETRY  INFORMATION  FOR  A(l, 1,2,1) 

AND  A ( 1, 1,2,2)  ETC. 

IF  ITA  =  1  TILT  SYMMETRY  IS  EVEN  WITH  RESPECT  TO  Z  SYM 

IF  ITA  =  2  TILT  SYMMETRY  IS  ODD  WITH  RESPECT  TO  Z  SYM. 

ITB  SYMMETRY  POINTER  FOR  C  AND  D  ARRAYS 

ITC  SYMMETRY  POINTER  FOR  E  AND  F  ARRAYS 

X  X  COMPONENT  OF  POSITION 

Y  Y  COMPONENT  OF  POSITION 

Z  Z  COMPONENT  OF  POSITION 

Y2  Y**2 

Z2  Z*  *2 

R2  X* *2  +  Y*  *2  +  Z**2 
R  SQRT (R2 ) 

I  DO  LOOP  VARIABLE.  IN  THE  FIELD  EXPANSION  LOOP  IT 
REPRESENTS  THE  POWER  TO  WHICH  X  IS  CURRENTLY  RAISED 
I.E.  X**(I-1) 

J  DO  LOOP  VARIABLE.  ALSO  Y**(2*J-2) 

K  DO  LOOP  VARIABLE.  ALSO  Z**(K-1) 

XB  X**(I-1) 

YEXB  X** (1-1) *Y** (2*J-2) 

ZEYEXB  X** (1-1) *Y* * (2 *J-2 ) *Z**  <K-1) 

IJK  I  +  2*J  +  K 

II  POINTS  TO  THE  ARRAY  LOCATION  WHERE  THE  CURRENT  POWER 
SERIES  COEFFICIENT  FOR  BX  IS  LOCATED 

JJ  BY  COEFFICIENT  LOCATION  POINTER 

KK  BZ  COEFFICIENT  LOCATION  POINTER 

BX, BY,BZ  ARE  USED  TO  CONSTRUCT  THE  MAGNETIC  FIELD  WITHIN  THE 
POWER  SERIES  LOOP. 

EXPR  EXP (- . 06*R2) 

TILTL  HOLDS  THE  LAST  VALUE  OF  THE  TILT  FOR  WHICH  THE  TILT 
INDEPENDENT  COEFFICIENTS  A-F  WERE  CALCULATED 
TT  A  REAL  ARRAY  HOLDING  THE  POWERS  OF  THE  TILT. 

TT ( 1) =TILT*  *0,  TT (2) =TILT**1,  ETC. 

CON  =0  FOR  R  LESS  THAN  2 

=1  FOR  R  GREATER  THAN  2.5 

GOES  FROM  0  TO  1  IN  THE  REGION  2  TO  2.5 

FOR  MORE  INFORMATION  CALL  OR  WRITE  K.  A.  PFITZER  OR  W .  P.  OLSON 
AT  MCDONNEL  DOUGLAS  ASTRONAUTICS  CO.  5301  BOLSA  AVE,  HUNTINGTON 
CALIF.,  PHONE  (714)  896-3231. 
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DIMENSION  BF(3) ,XX(3) , AA(64) , BB(64) , CC(44) , DD(44) ,EE(64) ,FF(64) , 

A (32) , B { 32 ) ,C(22) , D (22 ) ,E(32) , F ( 32 ) ,TT(4) , ITA(32) , ITB (22) , ITC (32) 

DATA  (ITA(I) ,1=1, 32)  /2, 1,2, 1,2, 2, 1,2, 1,2, 1,2, 1,2, 1,2, 2, 1,2, 2, 2,1, 

2, 1,2, 1,2, 1,2, 2, 2,1/ 

DATA  (ITB (I), 1=1, 22)  /2, 1,2, 1,2, 2, 1,2, 2, 2, 1,2, 1,2, 1,2, 1,2, 2, 2, 1,2/ 

DATA  (ITC  (I) ,1  =  1, 32)  /l, 2, 1,2, 1,1, 2, 1,2, 1,2,1, 2, 1,2,1, 1,2, 1,1, 1,2, 

1212121112/ 

DATA  '(AA  (I)  I  =  l!  64)  /-2.26836E-02,  -1.01863E-04,  3.42986E+00,  TOTAL 

-3. 12195E-04,  9 . 50629E-03, -2 . 91512E-06, -1 . 57317E-03,  8.62856E-08,  TOTAL 
-4 .26478E-05,  1 . 62 92 4E-08 , -1 . 2754 9E-04,  1 . 90732E-06, -1 . 65983E-02,  TOTAL 
8 . 46680E-09, -5 . 55850E-05,  1.37404E-08,  9.91815E-05,  1.59296E-08,  TOTAL 
4 . 52864E-07, -7 . 17669E-09,  4.98627E-05,  3 . 33662E-10, -5 . 97620E-02,  TOTAL 
1 . 60669E-05, -2 .294S7E-01,-1. 43777E-04,  1 . 09403E-03, -9 . 15606E-07,  TOTAL 
1 .60658E-03,-4.01198E-07,-3.15064E-06,  2.03125E-09,  4.92887E-04,  TOTAL 
-1 .80676E-07, -1 . 12022E-03,  5 . 98568E-07, -5 . 90009E-06,  5.16504E-09,  TOTAL 
-1 . 48737E-06,  4 . 83477E-10, -7 . 44379E-04,  3.82472E-06,  7.41737E-04,  TOTAL 
-1 . 31468E -05,-1. 2472 9E-04,  1. 92930E-08, -1 . 91764E-04, -5 . 30371E-08,  TOTAL 
1 . 38186E-05, -2 . 81594E-08,  7.46386E-06,  2.64404E-08,  2.45049E-04,  TOTAL 
-1 . 81802E-07, -1 . 00278E-03,  1 . 98742E-06, -1 . 1 6425E-05,  1.17556E-08,  TOTAL 
-2 . 46079E-06, -3.45831E-10,  1. 02440E-05, -1 . 9071 6E-08, -4 . 00855E-05,  TOTAL 
1 .25818E-07/ 

DATA  (BBC)  ,  1  =  1,  64)  /  9.47753E-02,  1 . 45981E-04 ,  - 1 . 82  93  3E+00 ,  TOTAL 

5 . 54882E-04,  5 . 03665E-03, -2 . 07698E-06,  1 . 10959E-01, -3 . 45837E-05,  TOTAL 
-4 . 40075E-05,  5 . 06464E-07, -1 .20112E-03,  3.64911E-06,  1.49849E-01,  TOTAL 
-7 . 44929E-C5,  2.46382E-04,  9 . 65870E-07, -9 . 54881E-04,  2.43647E-07,  TOTAL 
3 . 0652OE-34,  3.07836E-07,  6.48301E-03,  1 .26251E-06, -7 . 09548E-03,  TOTAL 
-1 . 55596E-05,  3 . 06465E+00, -7 . 84893E-05,  4.95145E-03,  3.71921E-06,  TOTAL 
-1 . 52002E-01,  6. 81988E-06, -8 . 55686E-05, -9 . 01230E-08, -3 . 7 1458E-04 ,  TOTAL 
1 . 30476E-07, -1 . 92971E-01,  1 . 51390E-05, -1 . 45912E-04 , -2 . 22778E-07 ,  TOTAL 

6.49278E-05,-3.72758E-08,-1.59932E-03,  8.04921E-06,  5.38012E-01,  TOTAL 
-1 . 43182E-04,  1  .  50000E-04,  5 . 88020E-07, -1 . 59000E-02,  1.60744E-06,  TOTAL 
3 . 17837E-04,  1 . 7 8959E-07 , -8 . 937 94E-03,  6.37549E-06,  1.27887E-03,  TOTAL 
-2 .45878E-07, -1 . 93210E-01,  6. 91233E-06, -2 . 80637E-04, -2 . 57073E-07,  TOTAL 
5 . 78343E-05,  4.52128E-10,  1 . 8 962  IE-04 , -4 . 84 91 IE-08, -1 . 50058E-02 ,  TOTAL 
6.21772E-C6/ 

DATA  (CC(I) ,  1  =  1,  44) /- 1 . 88 177E-02 , - 1 . 924 93E-0 6, -2 . 8  90  64E-01 ,  TOTAL 

-8 . 49439E-05, -4 . 76380E-04, -4 . 52998E-08,  1.61086E-03,  3.18728E-07,  TOTAL 
1.29159E-06,  5 . 52259E-10,  3.95543E-05,  5.61209E-08,  1.38287E-03,  TOTAL 
5 . 7  42  37E-07,  1.86489E-06,  7.10175E-10,  1 . 45243E-07, -2 . 97591E-10,  TOTAL 
-2 . 43029E-03, -6 . 70000E-07, -2 . 30624E-02, -6 . 22193E-06, -2 . 40815E-05,  TOTAL 
2 . 01689E-08,  1.76721E-04,  3.78689E-08,  9.88496E-06,  7.33820E-09,  TOTAL 
7 . 32126E-05,  8 . 43986E-08,  8 . 82449E-06, -6 . 1 1708E-08,  1.78881E-04,  TOTAL 
8 . 62589E-C7,  3.43724E-06,  2 . 53783E-09, -2 . 04239E-07,  8.16641E-10,  TOTAL 
1 . 6  8  0  7  5  E - 0  5 ,  7.62815E-09,  2.26026E-04,  3.66341E-08,  3.44637E-07,  TOTAL 
2.25531E-10/ 

DATA  (DD  (I) ,  1  =  1  ,  44) /  2.50143E-03,  1.01200E-06,  3.23821E+00,  TOTAL 

1 . 0858  9E-C5, -3 . 391 99E-05 , -5 . 27052E-07 , -  9 . 4 61 61E-02 , -1 . 95413E-09,  TOTAL 
-4 .23614E-06,  1 . 43153E-08, -2 . 62948E-L4,  1 . 05138E-07, -2 . 15784E-01,  TOTAL 

-2 .20717E-07, -2 . 65687E-05,  1.2637CE-08,  5 . 88917E-07, -1 . 13658E-08,  TOTAL 
1 . 64  385E-0  3 ,  1 . 44263E-06, -1 . 66045E-01, -1 . 46096E-05,  1.22811E-04,  TOTAL 
3 . 43922E-08,  9 . 66760E-05, -6 . 32150E-07, -4 . 97400E-05, -2 .78578E-08,  TOTAL 
1 .77366E-02,  2. 05 4 01E- 07,-1. 91756E-03, -9 . 4 93 92E-07 , - 1 . 99488E-01,  TOTAL 
-2 .07170E-06, -5 . 40443E-05,  1.59289E-08,  7.30914E-05,  3.38786E-08,  TOTAL 
-1 . 59537E-04, -1 . 65504E-07,  1.90940E-02,  2.03238E-06,  1.01148E-04,  TOTAL 

5 .20815E-08/ 

DATA  (EE(I) , I  C, 64) /-2 . 77 92 4E+0 1 , -1 . 01457E-03,  9.21436E-02,  TOTAL 

-8 . 52 177E-06,  5.19106E-01,  8 . 28881 E-05, -5 . 5965  IE-04 ,  1.16736E-07,  TOTAL 
-2 . 1 1206E-03, -5 . 35469E-07,  4 . 4 1 990E-0 1 , - 1 . 3367 9E-05, -7 . 18642E-04,  TOTAL 
6 . 17358E-08, -3.51 990E-03, -5 .29070E-07,  1 . 88443E-06, -6 . 60696E-10,  TOTAL 
-1 . 34708E-03,  1.02160E-07,  1.58219E-06,  2.05040E-10,  1.18039E+00,  TOTAL 
1 . 58903E-0  4 ,  1 . 86944E-02, -4 . 46477E-06,  5.49869E-02,  4.94690E-06,  TOTAL 
-1 . 18335E-04,  6 . 95684E-09, -2 .73839E-04, -9. 17883E-08,  2.79126E-02,  TOTAL 
-1 . 02567E-0'-  1.25427E-04,  3 . 07 1  43E-08, -5 . 31 82 6E-04  , -2 . 9847  6E-08,  TOTAL 
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*-4 . 89899E-05,  4.91480E-08,  3.85563E-01,  4.16966E-05,  6.74744E-04, 
* -2 . 087  3  6E-07,  -3 . 42  654E-03, -3 . 13  957E-06, -6 . 31361E-06, -2 . 92  981E-09, 
* -2 . 63883E-03, -1 . 32235E-07, -6 . 1940  6E-0  6,  3 . 5  4  334E-09,  6 . 65986E-03, 
*-5 . 81949E-06, -1 . 88809E-04,  3 . 62055E-08, -4 . 64380E-04, -2 .21159E-07, 
*-l . 77  4  96E-04 ,  4 . 95560E-08, -3 . 18867E-04, -3 . 17697E-07, -1 . 05815E-05, 

*  2 . 22220E-09/ 

DATA  (FF(I) , 1=1,64) /-5.07092E+00,  4 . 71960E-03, -3 . 79851E-03, 

*-3 . 67309E-06, -6 . 0243 9E -01,  1.08490E-04,  5.09287E-04,  5.62210E-07, 

*  7 . 05718E-02,  5 . 13160E-06, -2 . 85571E+00, -4 . 31728E-05,  1.03185E-03, 

*  1 . 05332E-07 ,  1.04106E-02,  1.60749E-05,  4.18031E-05,  3.32759E-08, 

*  1 . 20113E-01,  1.40486E-05,-3.37993E-05,  5.48340E-09,  9.10815E-02, 
* -4 .00608E-04,  3 .75393E-03, -4 . 69939E-07, -2 . 48561E-02,  1 . 31836E-04, 
*-2 . 67755E-04, -7 . 60285E-08,  3.04443E-03, -3 .28956E-06,  5.82367E-01, 

*  5.39496E-06, -6. 15261E-04,  4.05316E-08,  1 . 13546E-02, -4 . 26493E-06, 
*-2 . 72007E-02,  5 . 72523E-08, -2 . 98576E+00,  3.07325E-05,  1.51645E-03, 

*  1 . 25098E-06,  4.07213E-02,  1.05964E-05,  1.04232E-04,  1.77381E-08, 

*  1.92781E-01,  2 . 15734E-05, -1 . 65741E-05, -1 . 8868 3E -09,  2.44803E-01, 

*  1.51316E-05,-3.01157E-04,  8.47006E-08,  1 . 86971E-02, -6 . 94074E-06, 

*  9 . 13198E-03, -2 . 38052E-07,  1.28552E-01,  6 . 92595E-06, -8 . 36792E-05, 
*-6 . 10021E-08/ 

DATA  TILTL/99 . / 

SET  UP  SOME  OF  THE  INITIAL  POSITION  VARIABLES 
X=XX(1) 

y=xx<2) 

Z=XX (3) 

Y2=Y*  *2 
Z2=Z**2 
R2=X**2+Y2+Z2 

SET  MAGNETIC  FIELD  VARIABLES  TO  ZERO 
BX=0. 

BY=0  . 

BZ  =  0  . 


CHECK  TO  SEE  IF  POSITION  IS  WITHIN  REGION  OF  VALIDITY 
CON=l . 

IF  DISTANCE  TOO  LARGE  TAKE  ERROR  EXIT 
IF (R2 .GT.225 . )  GO  TO  50 

IF  DISTANCE  TOO  SMALL  SET  FIELD  TO  ZERO  AND  EXIT 

IF (R2 . LT . 4 . ) GO  TO  40 

IF  (R2 .LT. 6.25)  CON=CON* (R2-4 . 0) /2 . 25 

IF  TILT  HAS  NOT  CHANGED,  GO  DRECTLY  TO  FIELD  CALCULATION 

IF (TILTL . EQ . TILT) GO  TO  6 

SET  UP  POWERS  OF  TILT 

TILTL=TILT 

TT (1) =1 

TT (2) =TILTL 

TT (3) =TILTL**2 

TT (4) =TILT*TT  (3) 

SET  UP  THE  X  AND  Z  COMPONENT  TILT  INDEPENDENT  COEFFICIENTS 
DO  1  1=1,32 
J= ( I - 1 ) *2+1 
K=ITA ( I ) 

A { I )  =AA ( J) *TT (K) +AA ( J+l )  *7T(K  +  2) 

B ( I ) =BB ( J) *TT (K) +BB (J+l ) *TT(K+2) 

K= ITC ( I ) 

E ( I )  =EE ( J) *TT (K) +EE ( J  +  l ) *TT(K  +  2) 

F  (I)  =FF  (J)  *TT  (K)  +FF  (J  +  l)  *TT  (K  +  2) 

1  CONTINUE 

SET  UP  THE  Y  COMPONENT  TILT  INDEPENDENT  COEFFICIENTS 


TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 
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DO  2  1=1,22 
J=(I-1) *2+1 
K=ITB (I) 

C  (I)  =CC  ( J)  *TT (K) +CC (J  +  l) *TT (K+2) 

D (I) =DD  (J)  *TT  (K)  +DD  (J+l)  *TT  (K  +  2) 

2  CONTINUE 

6  £XPR=EXP (-0 . 06* R2) 

C 

INITIALIZE  THE  POINTERS 

11  =  1 

JJ=1 

KK=1 

XB=1  . 

BEGIN  SUM  OVER  X 
DO  30  1=1,5 
YEXB=XB 

BEGIN  SUM  OVER  Y 
DO  20  J=l, 3 

IF ( I +2 * J . GT .  8)  GO  TO  25 
ZEYEXB=YEX3 
I JK=I+2  * J+l 
K  =  1 

Z  LOOP  STARTS  HERE 

10  BZ=BZ+ (E (KK) +F (KK) *EXPR) *ZEYEXB 
KK=KK+1 

8X=BX+ ( A ( I I ) +B (II) *EXPR  ) *ZEYEXB 
11=11+1 

IF ( I JK  .GT.  8)  GO  TO  15 

BY=BY+ (C ( JJ) +D ( JJ) *EXPR  ) *ZEYEXB*Y 

JJ=JJ+1 

ZEYEXB=ZEYEXB*Z 
I JK=IJK+1 
K=K+ 1 

IF ( I JK . LE . 9 . AND . K . LE . 5 )  GO  TO  10 
15  YEXB=YEXB*Y2 
20  CONTINUE 
25  XB=XB*X 
30  CONTINUE 

SET  UP  THE  OUTPUT  ARRAY,  MULTIPLY  BY  CON.  CON  IS  NORMALY  ONE 
BUT  INSIDE  OF  R=2 . 5  IT  GOES  TO  ZERO.  INSIDE  R=2  IT  IS  ZERO. 

4  0  BF ( 1 ) =BX  *CQM 

BF (2) =BY*CON 
BF (3) =BZ*CON 
RETURN 

ERROR  EXIT  IF  OUTSIDE  OF  R  =  15. 

50  WRITE (6, 60)  XX 

00  FORMAT  ( 4H  X=  ,E10.3,4H  Y=  ,E10.3,4H  Z=  ,E10.3,’76H  IS  OUTSIDE  THE 
•VALID  REG I ON - -POWER  SERIES  DIVERGES  BFIELD  IS  SET  TO  ZERO  ) 

GO  TO  40 
END 
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Appendix  C 
Subroutine  INVARM 

This  appendix  presents  a  listing  of  subroutine  INVARM,  the  central  routine  for  developing 
the  new  radiation  belt  models.  The  operation  and  functions  performed  by  this  routine  are 
spelled  out  in  section  2.3.  The  routine  is  preceded  by  a  test  routine  that  presents  a  sample 
of  its  capabilities  and  also  provides  a  means  for  assessing  its  operation. 

The  test  routine  varies  the  distance,  latitude,  longitude,  and  universal  time  and  asks  the 
INVARM  program  to  calculate  the  invariant  for  18  different  pitch  angles  at  each  location. 

All  input  and  out  variables  are  passed  via  the  arguments  of  the  subroutine  call.  It  is 
possible  to  limit  the  number  of  terms  used  by  the  internal  field  routine  by  setting  NMAX  to  a 
value  between  2  and  10.  This  must  be  done  via  labeled  COMMON  /GCOM/.  The  internal 
field  routine  checks  the  value  of  this  variable.  If  it  is  set  to  zero  as  it  normally  is  when 
labeled  COMMON  is  preset  by  the  loader  or  set  to  any  value  other  than  2  through  10,  then 
the  internal  field  routine  uses  the  maximum  coefficients  defined  for  the  IGRF  field.  For  the 
coefficients  supplied  with  this  program  an  NMAX  of  1 1  is  used.  Labeled  COMMON 
/MOMENT/  contains  the  dipole  moment  of  the  main  field  as  calculated  from  the  coefficient 
set  that  is  in  use.  This  may  be  used  by  any  routine  that  requires  it. 

The  calling  argument  for  the  routine  are 

Input  variables 

XLAT  The  geographic  latitude  measured  in  degrees  from  the  geographic  equator, 
plus  is  north  and  minus  is  south. 

XLONG  The  geographic  east  longitude  measured  from  Greenwich  England  (0-360) 

R  The  radial  distance  from  the  center  of  the  earth  in  unit  of  earth  radii.  One 

radius  is  6371 .2  km. 

YR  The  year  of  the  calculation.  This  variable  is  used  by  the  internal  magnetic 

field  routine  to  calculate  the  epoch  of  the  magnetic  field.  Whenever  this  value 
changes  by  0.1  years  the  coefficient  set  for  the  internal  field  is  updated.  It  is 
suggested  that  this  value  not  be  changed  unless  the  drift  of  the  internal  field 
during  the  calculations  is  important  to  the  analysis.  The  coefficients  are  valid 
from  1945  to  the  present.  Dates  earlier  than  1945,  will  cause  the  routine  to 
use  the  1945  coefficient  set.  Use  caution  in  predicting  the  field  far  into  the 
future.  Historically,  predicting  the  field  into  the  future  has  not  been  very 
successful. 

DAY  The  day  of  the  year.  January  1  is  DAY  =  1 .  This  is  a  floating  point  variable 

but  it  should  be  limited  to  whole  numbers. 
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TIME  The  universal  time  in  hours.  This  is  a  floating  point  number  should  represent 

the  correct  universal  time  to  the  required  precision. 

JSWITCH  An  integer  variable  that  controls  whether  the  external  field  is  included  in  the 
calculation.  JSWITCH  negative  uses  the  internal  field  only,  JSWITCH  =  0  or 
positive  uses  the  internal  plus  external  magnetic  field. 

NUMANG  An  integer  variable  that  specifies  the  number  of  pitch  angles  for  which  the 
invariant  must  be  calculated. 

PANGLE  A  single  variable  or  an  array  that  contains  the  pitch  angles  for  which  the 
routine  is  to  calculate  the  invariant.  The  dimension  of  PANGLE  must  be 
equal  to  or  greater  than  NUMANG.  If  NUMANG  is  1 ,  then  PANGLE  may  be  a 
simple  undimensioned  variable. 

Output  parameters 

EL  A  simple  variable  or  an  array  dimensioned  to  at  least  NUMANG  that  will 

contain  the  L  value  for  the  specified  location  and  pitch  angle.  If  no  L  could  be 
calculated  EL  is  set  to  -1 .0,  if  the  mirror  point  is  below  1 .03  Re  or  El  is  set  to 
1 00  if  the  field  line  is  open  or  the  maximum  number  of  steps  is  reached  by 
the  routine. 

BLOCAL  The  value  of  the  magnetic  field  at  the  observation  point. 

BMIN  The  minimum  value  of  the  magnetic  field  along  the  particle  line  of  force. 

XMLONG  The  magnetic  longitude  of  the  minimum  B  point  on  the  magnetic  line  of  force. 

0  degrees  is  local  midnight. 

XMLAT  The  magnetic  latitude  of  the  observation  point. 

BMAXAN  The  Mirror  point  magnetic  field  for  each  of  the  pitch  angles.  BMAXAN  can 
either  be  a  simple  variable  or  and  array  dimensioned  at  least  to  order 
NUMANG. 

XJ  The  value  of  the  second  invariant  for  each  pitch  angle.  XJ  can  either  be  a 

simple  variable  or  and  array  dimensioned  at  least  to  order  NUMANG. 

DENSTY  The  column  density  of  the  atmosphere  along  the  particle's  bounce  path  in 

gm/cm2.  DENSTY  can  either  be  a  simple  variable  or  and  array  dimensioned 
at  least  to  order  NUMANG. 
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COMMON/GCOM/  ST, CT, SP,  CP,  AOR, BT, BP, BR, NMAX, YEARI 

DIMENSION  BMAXAN (20) ,XJ (20) , ANGLE (20) ,EL(20) ,DENSTY(20) ,ALPHEQ(20) 
common/temp/nlast, n21ast 
C 

C  THIS  PROGRAM  PROVIDES  A  TEST  RUN  OP  THE  L  VALUE  SUBROUTINES 

C 

CHARACTER* 6  IAR(3) 

DATA  IAR/6H  INT  , 6H  IN+EX, 6H  L  AVE/ 

DO  500  1=1, 18 
ANGLE (I) =90- < I- 1 ) *5 

500  CONTINUE 
NUMANG= 1 8 
IDSWIT=1 
LN=100 
YEAR=1990 
DA=1 

DO  5  IU=1, 1, 12 
UT=IU-1 
LN=1 00 

DO  4  IL=1, 31, 30 
FLAT=IL-1 

DO  3  ILG=1, 181, 180 

XLONG=ILG 

DO  2  IR=2, 8, 2 

R=IR- . 5 

DO  1  IC-1,2 

call  gettim(ihr, imin, isec, 1100) 
btime=float (imin*60+isec) +float (ilOO)  /100 . 

CALL  INVARM (FLAT, XLONG, R, YEAR, DA, UT, IC-2 , ANGLE, NUMANG, 

*EL, BLOC, BM, XMLONG, XMLAT, BMAXAN, XJ, DENSTY) 
call  gettim(jhr, jmin, jsec, jlOO) 
tine= float ( jmin* 60+ jsec) +f loat (jl00)/100. -btime 
IF ( IC . EQ . 1) WRITE ( * , 103) 

103  FORMAT (1H1) 

WRITE (*,  101) FLAT,  XLONG, R, YEAR, DA, UT, IAR(IC) , BLOC, BM, XMLAT, XMLONG 

101  FORMAT  (II,'  Lat  =  ',f6.1,'  Long  =  ',f7.1,'  R  =  \f4.1, 

*/,  '  Year  =  \f7.1,'  Day  =  ’,f5.0,’  UT  =  ',f6.2,'  Field  =',A6, 
*/,  '  Blocal  =  ' , F8 . 5, '  Bmin  =  ',F8.5,'  Mlat  =  ',f8.3, 

*'  Mlong  =  ', f9.3) 

WRITE (*, 102) 

102  FORMAT(/,'  P.  Angle  B  mir  2nd  Inv.  L  Density', 

*'  Eq.  Pitch  Angle',/) 

LN=0 

DO  50  1=1, NUMANG 

50  ALPHEQ ( I ) =ASIN (SQRT (BM/BLOC) *SIN ( ANGLE ( I ) *.01745329))/. 0174532 9 

write  ( *, 100)  (angle ( i) , bmaxan (i),xj(i),el(i), densty (i) , 

*alpheq(i) , i=l,numang) 

10  wr ite ( * , *) nlast, n21ast , time 

100  format (0Pf6 . 1 , 2f 11 . 5, f8 . 3, 1PE15 . 5, 0PF13 . 2 ) 

1  CONTINUE 

2  CONTINUE 

3  CONTINUE 

4  CONTINUE 

5  CONTINUE 
END 
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SUBROUTINE  INVARM (XLAT, XLONG, R, YR, DAY, TIME, JSWTCH, PANGLE, 

•NUMANG, EL, BLOCAL, BMIN, XMLONG, XMLAT, 8MAXAN, XJ, DENSTY) 

PURPOSE 

CALCULATE  THE  VARIOUS  MAGNETIC  COORDINATES  OF  THE  PARTICLES' 

DRIFT  SHELL.  CALCULATE  THE  1ST  AND  2ND  ADIABATIC  INVARIANTS  AND 
THE  L  PARAMETER  FOR  A  NUMBER  OF  PITCH  ANGLES  AT  THE  SPECIFIED 
LOCATION.  ALSO  DETERMINE  THE  LOCAL  MAGNETIC  FIELD,  THE  MAGNETIC 
LATITUDE,  THE  MINIMUM  MAGNETIC  FIELD  ON  THE  FIELD  LINE  AND  THE 
MAGNETIC  LONGITUDE  AT  THE  FIELD  MINIMUM. 

INPUT  —  ARGUMENT  LIST 

XLAT  GEOCENTRIC  GEOGRAPHIC  LATITUDE  IN  DEGREES  (+  IS  NORTH) 
XLONG  GEOCENTRIC  GEOGRAPHIC  LONGITUDE  EAST  OF  GREENWHICH  IN 
DEGREES 

R  GEOCENTRIC  DISTANCE  FROM  THE  EARTHS  CENTER  IN  UNITS 

EARTH  RADII,  RE.  RE=6371.2  KM 

YR  THE  YEAR  -  USED  BY  THE  INTERNAL  MAGNETIC  FIELD  ROUTINE 

TO  TAKE  INTO  ACCOUNT  THE  SECULAR  VARIATIONS 
(E.G.  JULY  15,1964  =  1964.54) 

NOTE***  YR  SHOULD  BE  CHANGED  ONLY  EVERY  FEW  DAYS  OR 
MONTHS.  NEW  FIELD  COEFFICIENTS  MUST  BE  COMPUTED  FOR 
EVERY  CHANGE  IN  YR,  THIS  COULD  CAUSE  A  LARGE  INCREASE  IN 
COMPUTER  TIME.  THE  EARTHS  FIELD  CHANGES  ONLY  ABOUT 
.001  GAUSS/YEAR  AT  THE  EARTHS  SURFACE. 

IF  YR  IS  CHANGED  BY  MORE  THAN  .1  YEAR  NEW  FIELD  COEFFS . 

ARE  COMPUTED 

DAY  THE  DAY  OF  YEAR  (1.-366.).  THE  DAY  IS  USED  BY  THE 
MAGNETIC  FIELD  ROUTINE  TO  CALCULATE  THE  TILT  OF  THE 
DIPOLE  AXIS  FOR  THE  EXTERNAL  FIELD  ROUTINE 
DAY  MUST  BE  A  WHOLE  NUMBER  AND  DAY  1  IS  JANUARY  1 
TIME  UNIVERSAL  TIME  IN  HOURS  (0.000-24.0000) 

JSWTCH  A  FLOW  CONTROL  VARIABLE 

JSWTCH  =-l  OR  NAEGATIVE,  COMPUTE  L  USING  INTERNAL 
FIELD  ONLY 

=  0  OR  POSITIVE,  COMPUTE  L  USING  INTERNAL  PLUS 
EXTERNAL  FIELD 

PANGLE  A  SINGLE  PITCH  ANGLE  OR  AN  ARRAY  OF  PITCH  ANGLES  FOR  THE 
INVARIANTS  WILL  BE  CALCULATED.  THE  PITCH  ANGLE  MUST  BE 
. LE.  90  AND  GT.O  AND  THE  ARRAY  MUST  BE  ORDERED  IN 
DESCENDING  ORDER  (90,  80,  70,...) 

NUMANG  THE  NUMBER  OF  ELEMENT  IN  THE  PANGLE  ARRAY 
OUTPUT  PARAMETERS 

EL  A  SINGLE  VARIABLE  OR  AN  ARRAY  OF  DIMENSION  NUMANG.  THIS 

RETURNS  THE  L  VALUE  CALCLUATED  FROM  THE  INVARIANT. 

VARIABLE  JSWTCH 
*****  NOTE  *  *  *  * 

SINCE  THIS  ROUTINE  USES  AN  ACTUAL  MAGNETOSPHERIC 
MAGNETIC  FIELD,  THE  FIELD  LINES  ARE  NOT  ALL  CLOSED. 

THUS  L  IS  DEFINED  ONLY  IN  THE  INNER  MAGNETOSPHERE  (IN 

THE  REGION  OF  CLOSED  DRIFT  SHELLS) .  AN  ATTEMPT 

TO  CALCULATE  L  OUTSIDE  OF  THIS  REGION  WILL  SET  EL  TO 

100  (EL- 100),  SET  BMIN  TO  THE  LOCAL  FIELD  VALUE 

AND  SET  XMLONG  TO  ZERO  UNLESS  MINIMUM  B  WAS  PASSED  PRIOR 

TO  THE  DETECTION  OF  THE  ERROR 

IF  THE  MIRROR  POINT  FOR  A  GIVEN  PITCH  ANGLE  IS  BELOW  200KM 
THEN  EL  IS  SET  TO  -1 

BLOCAL  THE  VALUE  OF  THE  MAGNETIC  FIELD  AT  THE  INPUT  POSITION 

(IN  GAUSS) 

BMIN  THE  MINIMUM  VALUE  OF  B  ALONG  THE  FIELD  LINE  IN  GAUSS 
XMLONG  THE  MAGNETIC  LONGITUDE  OF  THE  MAGNETIC  FIELD  MINIMUM 
MEASURED  EAST  OF  THE  PRIME  MAGNETIC  MERIDIAN 

( IN  DEGREES) . 

A  FT  CONSTANT  IN  SUBROUTINE  MG  LONG  ALLOWS  THE  USER 
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TO  SELECT  EITHER  A  CENTERED  DIPOLE  MAGNETIC  COORDINATE 
SYSTEM  WITH  ZERO  AT  69  DEG  W.  GEOGRAPHIC,  OR  AN  OFFSET 
DIPOLE  COORDINATE  SYSTEM  WITH  ZERO  THROUGH  GREENWHICH. 
XMLAT  THE  MAGNETIC  LATITUDE  IN  DEGREES  OF  THE  CURRENT  POSITION 
BMAXAN  A  SINGLE  VARIABLE  OR  AN  ARRAY  OF  AT  LEAST  DIMENSION  NUMANG 
THAT  WILL  HOLD  THE  MIRROR  POINT  MAGNETIC  FIELD  FOR  THE 
VARIOUS  PITCH  ANGLE. 

XJ  A  SINGLE  VARIABLE  OR  AN  ARRAY  OF  AT  LEAST  DIMENSION  NUMANG 

THAT  WILL  HOLD  THE  VALUES  OF  THE  SECOND  INTEGRAL  INVARIANT 
FOR  EACH  PITCH  ANGLE 

CONSTANTS 

ERR  =  0.0005  SCALES  THE  ERROR  LIMITS  FOR  THE  INTEGRATION 
THE  ERROR  IN  L  IS  APPROXIMATELY  L*ERR 

SUBROUTINES  REQUIRED 
SUBROUTINE  STEPSZ 
SUBROUTINE  BMNEXT 
SUBROUTINE  HILTEL 
SUBROUTINE  INTERP 
SUBROUTINE  INTGRT 
SUBROUTINE  INVR 
SUBROUTINE  MGLONG 

VARIABLES  (PARTIAL  LIST  TO  HELP  UNDERSTAND  THE  CODE) 

BINTL  A  REAL  ARRAY  THAT  SAVES  THE  VALUE  OF  THE  MAGNETIC  FIELD 
AT  THE  INPUT  POSITION 

B  A  REAL  ARRAY  THAT  HOLDS  THE  INSTANTANEOUS  MAGNETIC  FIELD 

VECTOR  AT  EACH  INTEGRATION  STEP 
BB  A  2  DIMENSIONED  REAL  ARRAY  THAT  HOLDS  ALL  OF  THE 

MAGNETIC  FIELD  MAGNITUDES  CALCULATED  AT  ALL  OF  THE 
INTEGRATION  STEPS 

BL  A  REAL  ARRAY  THAT  SAVES  THE  MAGNETIC  FIELD  VECTOR 

FROM  THE  PREVIOUS  INTEGRATION  STEP 
BMAX  THE  MAGNETIC  FIELD  AT  THE  PARTICLE  MIRROR  POINT 
COSINE  THE  COSINE  OF  THE  GEOGRAPHIC  LATITUDE 
DAYYR  THE  DAY  OF  THE  YEAR 

DDS  THE  ESTIMATED  STEP  SIZE  NECESSARY  TO  COMPLETE  THE 
INTEGRATION.  IF  NO  ESTIMATE  IS  YET  POSSIBLE  IT  IS 
SET  TO  100. 

DEL  SCALES  THE  INTEGRATION  STEP  SIZE.  IT  IS  PROPORTIONAL 
TO  THE  FOURTH  ROOT  OF  THE  ERROR  LIMITS.  IF  IT  IS 
POSITIVE  INTEGRATION  WILL  BE  PARALLEL  TO  THE  FIELD,  IF 
NEGATIVE  IT  IS  ANTIPARALLEL 

DS  THE  CURRENT  VALUE  OF  THE  INTEGRATION  STEP  SIZE  IN 

EARTH  RADII.  POSITIVE  IS  FOR  PARALLEL  TO  FIELD, 

NEGATIVE  FOR  ANTIPARALLEL 

IERFLG  AN  ERROR  FLAG  SET  BY  SUBROUTINE  INTGRT.  IF  NON-ZERO 
THE  INTEGRATION  HAS  GONE  BEYOND  THF,  SET  LIMITS  AND 
MUST  BE  TERMINATED 

JSW  A  FLOW  CONTROL  PARAMETER  USED  BY  THE  MA.GNETIC  FIELD 
SUBROUTINE 

KODE  SET  EQUAL  TO  ONE  TO  INDICATE  TO  SUBROUTINE  BMNEXT 
THAT  CARESIAN  COORDINATES  ARE  TO  BE  USED 
KS  A  VARIABLE  TRANSMITTED  TO  SUBROUTINE  INTERP.  IT  IS 

USED  TO  DETERMINE  WHICH  SOLUTION  APPLIES  TO  THE 
PARTICULAR  INTERPOLATION 

MINFLG  INITIALLY  SET  TO  ZERO.  IT  IS  SET  TO  ONE  WHEN  THE  FIELD 
MINIMUM  HAS  BEEN  PASSED 
N  THE  CURRENT  INTEGRATION  STEP  NUMBER 

PICON  PI  /  180. 

Q,  QL  REAL  ARRAYS  CONTAINING  THE  CURRENT  AND  PREVIOUS  ERROR 

ESTIMATES.  USED  BY  GILLS  METHOD  INTEGRATION  ROUTINE  TO 
CONTROL  ROUND  OFF  ERROR 
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RMIN  THE  VECTOR  POSITION  TO  THE  MAGNFTIC  MINNIMUM 

RMAG  THE  MAGNITUDE  OF  THE  DISTANCE  TO  BMIN 

SER  ERROR  CONTROL  VARIABLE.  THE  INTEGRATION  STOPS  IF 

THE  CURRENT  POSITION  POINT  IS  WITHIN  DISTANCE  SER  OF 
BMAX 

SF  OUTPUT  OF  THE  INTERPOLATION  SUBROUTINE  INTERP .  IT 

INDICATES  THE  SCALAR  DISTANCE  ALONG  THE  FIELD  WHERE 
B  IS  EQUAL  TO  BMAX 

SXJ  A  REAL  ARRAY  WHICH  SAVES  THE  INTEGRATION  STEP  VALUES 
OF  THE  SECOND  ADIABATIC  INVARIANT 
UT  UNIVERSAL  TIME 

XDS  TEMORARY  VALUE  USED  FOR  OBTAINING  DISTANCE  TO  COMPLETION 

OF  INTEGRATION 

XJ  FINAL  VALUE  OF  THE  SECOND  ADIABATIC  INVARIANT 

XL  A  REAL  ARRAY  HOLDING  THE  PREVIOUS  VALUE  OF  THE  POSITION 

VECTOR 

XSV  A  3  DIMENSIONED  REAL  ARRAY  HOLDING  ALL  OF  THE  POSITION 
VECTORS  ALONG  THE  INTEGRATION  PATH 
XXJ  INTERPOLATED  VALUE  OF  THE  SECOND  ADIABATIC  INVARIANT 
YEAR  THE  YEAR 

ZP  THE  Z  COMPONENT  OF  THE  POSITION  VECTOR  IN  CENTERED 

DIPOLE  COORDINATES 

VERSION  6/91 

FOR  MORE  INFORMATION  CALL  OR  WRITE  K.  A.  PFITZER  AT  MCDONNELL 
DOUGLAS  ASTRONAUTICS  CO.  5301  BOLSA  AVE,  HUNTINGTON  BEACH  CALIF. 
PHONE  (714)  896-3231. 

COMMON/BXYZCM/YEAR, DAYYR,  UT,  KODE,  JSW 

COMMON  /INTPAR/DS,DEL,N, IERFLG, XL (3) ,XSV(100, 3, 4) , 

*RSV (100) , RMIN (3) , RMAG, IDSW, 

*QL  <  3 ) , Q ( 3 ) , BL ( 3 ) , SXJ (100) ,DDS 
common/temp/nlast,  n21a^t 

DIMENSION  BB (100, 4) , BB2 (100, 4) ,B(3) ,B2 (3) ,X(3) ,X2 (3) ,S(100) , 

*S2 (100) , DEN (100) , DEN2  (100) 

DIMENSION  EL ( * ) , PANGLE ( * ) , BMAXAN ( * ) , XJ(‘) ,BLL(3) , BLL2 (3) 

DIMENSION  XX  (3) , BINTL (3) , DENSTY (*) 

DATA  PICON/ . 01745329252/ 

DATA  ERR/. 0005/ 

DATA  CONI/ . 95/ 

OBTAIN  THE  CARTESIAN  COMPONENTS  OF  THE  POSITION  VECTOR 


C  CHECK  THE  PITCH  ANGLES,  THEY  MUST  BE  BETWEEN  90  AND  0  AND  THEY  MUST 

C  BE  IN  DESCENDING  ORDER  (IE.  90,85,80,.... 

r» 

N  MAN  G = N  U  MAN  G 

CALL  CHECK ( F ANGLE, NMANG, IER) 

IF  (IER.GT.0)  THEN 

WRITE(*,*) 'Pitch  angle  error,  must  be  monotonic  and  >0  &  <=90' 

DO  5  1=1, NMANG 
XJ  ( I ) -- 1 
BMAXAN ( I ) = - 1 
EL ( I ) --1 
DENSTY ( I ) =- 1 
5  CONTINUE 

RETURN 
END  IF 

COS INE=COS (XLA?*F ICON) 

XX  ( 1)  =R*COSINF.*CG5  ( XLGNG ‘P ICON) 

XX (2) =R  *COS I NE  *S IN ( XLGNG* P ICON ) 

XX ( 3 ) =  R  *  S I N ( X  DAT  *  P I CON ) 

C 
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ROTATE  TO  DIPOLE  COORDINATES  (FIRST  ROTATE  ABOUT  Z  291  DEGREES 
THEN  ABOUT  THE  NEW  Y  11.7  DEGREES  TO  THE  DIPOLE  AXIS) 

ZP=(XX(1) * . 3583679495-XX (2) *.9335804265) * .2027872954 
*+XX (3) *.9792228106 

EVALUATE  THE  MAGNETIC  LATITUDE 
XMLAT«90.-ACOS(ZP/R) /PICON 

SET  THE  MAGNETIC  LONGITUDE  TO  ZERO.  IF  MINIMUM  B  IS  REACHED 
PRIOR  TO  AN  ERROR  BEING  DETECTED  XMLONG  IS  UPDATED  TO  REFLECT 
MAGNETIC  LONGITUDE  AT  MINIMUM  B 
XMLONG=0 . 

SET  UP  THE  COMMON  BLOCK  INPUT  VARIABLES  FOR  THE  MAGNETIC  FIELD 

SUBROUTINE 

YEAR=YR 

UT-TIME 

DAYYR-DAY 

JSW=JSWTCH 

KODE=l 

IBEFLG=0 

EVALUATE  THE  MAGNETIC  FIELD  AT  THE  STARTING  POINT 
CALL  BMNEXT ( XX , B , BB (2,1)) 

BLOCAL  =  BB (2,1) 

BB2 (2,1) =BB (2, 1) 

SAVE  THE  INITIAL  POSITION  AND  MAGNETIC  FIELD  VECTORS 
DO  10  1-1,3 
BINTL ( I ) =B ( I ) 

B2 ( I) -B ( I ) 

XL ( I ) =XX ( I ) 

X2 ( I ) =XX ( I ) 

XSV (2 , 1,1) =XX ( I ) 

10  CONTINUE 
RSV (2) =R 

EXIT  THE  ROUTINE  IF  POSITION  IS  OVER  THE  POLAR  CAP  OR  DISTANCE 
IS  TOO  LARGE  OR  MAGNETIC  FIELD  IS  TOO  WEAK 

IF (ABS (XMLAT) .GT.75. .OR.R.GT.12. ,OR.BB(2,l) . LT . .00025)  THEN 
NMANG=0 
BMAXAN ( 1 )  =  1 0  0  . 

GOTO  218 
END  IF 

SET  UP  THE  INITIAL  VALUES  FOR  THE  VARIABLES 
NLAST=2 
N2LAST=2 
S (2) =0 . 

S2 (2) =0  . 

DEN (2) =0 
DEN2 (2) =0 
DDS=1 00 . 

MINFLG=0 

SET  BMIN  TO  LOCAL  FIELD  VALUE.  IF  MINIMUM  B  IS  REACHED  PRIOR 
TO  ERROR  DETECTION  BMIN  IS  UPDATED  TO  MINIMUM  B. 

BMIN=BB ( 2 , 1) 

SET  UP  THE  ERROR  LIMITS  FOR  THE  INTEGRATION 
SER=SQRT (ERR) 

C  STEP  SIZE  GOES  AS  ERROR  TO  THE  .25  POWER 

DEL=-2 . 5  *ERR*  *.25 
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DS=SER 


STEP  ONCE  IN  THE  INCREASING  FIELD  DIRECTION  AND  SET  STEP 
PARAMETERS  TO  INTEGRATE  IN  THE  DECREASING  FIELD  DRECTION 
IFLAG=0 

I F ( XMLAT . GT . 0 . ) GO  TO  30 

DEL=-DEL 

DS=-DS 

N=2 

DO  31  1=1,3 
Q  ( I )  =0 
X ( I ) =XL ( I ) 

CONTINUE 

CALL  INTGRT  (X, B, BB, S, DEN) 

IFLAG=IFLAG+1 

IF ( (BB (3, 1 ) . LT .BB (2, 1) ) .AND. (IFLAG.LE. 1) )  GOTO  20 
S  ( 1 )  =  S  <  3 ) 

DEN (1) =DEN  (3) 

BB ( 1 , 1 ) =BB ( 3 , 1) 

DO  34  1=1,3 

XL { I ) =X ( I ) 

BL ( I ) =B ( I ) 

Q ( I ) =0 

X  ( I ) =XX ( I ) 

B (I) =BINTL(I) 

CONTINUE 

RSV ( 1 ) =SQRT ( XL ( 1 ) **2+XL(2) **2+XL<3) **3) 

DELSV=DEL 

BEGIN  THE  FIELD  LINE  INTEGRATION.  THE  INTEGRATION  USES  A  VARIABLE 
STEPSIZE  WHICH  IS  DEPENDENT  ON  THE  CURVATURE  OF  THE  FIELD  LINE 
AND  ON  THE  DISTANCE  EACH  POINT  IS  FROM  EARTH  CENTER  (A  MEASURE 
OF  THE  MAGNETIC  FIELD  STRENGTH) .  THE  INITIAL  INTEGRATION  IS  A 
LINE  INTEGRAL  OF  THE  MAGNETIC  FIELD  UNIT  VECTOR.  THIS  INTEGRATION 
LOOP  ALSO  SAVES  ALL  OF  THE  VARIABLES  WHICH  ARE  LATER  NEEDED  TO 
EVALUATE  THE  SECOND  INTEGRAL  INVARIANT. 

DO  216  IA= 1 , NMANG 

DEL=DELSV 

DDS=100 

BMAX=BB (2, 1) /SIN (PANGLE ( IA) *PICON) **2 
N=NLAST 

IF  ( IA.NE . 1 ) THEN 
DO  41  1=1, 3 
BL ( I ) =BLL { I ) 

CONTINUE 
DS=S <N) -S (N-l ) 

CALL  INTEPP (BB (N-2, 1) , S (N-2) , BMAX, SF, 3) 

IF (ABS (SF) .GT.AES (S (N) ) )  THEN 
XDS=SF-S ( N) 

IF(ABS(XDS) .LE.SER)GOTO  100 
DDS=CONI *XDS 
END  IF 
END  IF 

CALL  STEPSZ (X, B, BB) 

CALL  INTGRT (X, B, BB, 3, TEN) 

IF  FIELD  IS  STILL  DECREASING  RELOOP 
IF (BB (N, 1 ) . LT . BB (N-l , 1 ) )  GOTO  40 

IF  MINIMUM  VALUES  HAVE  BEEN  CALCULATED,  JUMP  OVER  MINIMUM  ROUTINES 
WHEN  THE  CURRENT  VALUE  OF  B  EXCEEDS  THE  LAST,  FIND  THE 
INTERPOLATE''  MINIMUM  MAGNETIC  FIELD  VALUE  AND  USE  THIS  VALUE  TO 
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UPDATE  THE  VALUE  OF  BMAX  TO  REFLECT  THE  AVERAGE  DRIFT  SHELL  (IF 
AVERAGE  SHELLS  ARE  REQUIRED) 

USE  THE  DISTANCE,  SF,  TO  THE  FIELD  MINIMUM  TO  DETERMINE  THE 
MAGNETIC  LONGITUDE  OF  THE  FIELD  MINIMUM 

IF (MINFLG.NE . 0)  GO  TO  50 
CALL  INTER? (BB(N-2, 1) ,S(N-2) , BMIN, SF, -1) 

CALL  MGLONG (XSV (N-2 , 1 , 1 ) , S (N-2 ) , SF, XMLONG, RMIN ( 1 ) , RMAG) 

MINFLG=1 

CONTINUE  STEPPING  ALONG  THE  FIELD  LINE  AS  LONG  AS  B  IS  LESS  THAN 
BMAX  AND  THE  INTEGRATION  IS  MORE  THAN  A  DISTANCE  SER  FROM  BMAX 
IF  BMAX  HAS  BEEN  EXCEEDED,  EXIT  TO  INTERPOLATION  SCHEME 

50  IF(BB(N, I) .GE.BMAX)  GO  TO  70 

IF  WE  ARE  OUTSIDE  OF  VALID  REGION  EXIT 
IF ( I ERF LG . NE . 0) THEN 
NMANG=IA-1 
IF ( IERFLG .GT . 0) THEN 
BMAXAN ( I A) =100 
ELSE 

BMAXAN ( IA) =-l 
ENDIF 
GOTO  218 
ENDIF 

CALL  INTERP  (BB (N-2, 1) , S (N-2) , BMAX, SF, 3) 

DDS=100 . 

IF  S  DOES  NOT  INCREASE  MONOTONICALLY,  IGNORE  INTERPOLATION 
AND  RELOOP 

IF (ABS ( SF) . LE . ABS ( S (N) ) )  GO  TO  40 
XDS=SF-S (N) 

IF  WITHIN  SER  OF  BMAX  STOP  INTEGRATION  GO  GET  VALUE  OF  INVARIANT 
IF(ABS(XDS) .LT.SER)  GO  TO  100 
DDS=CONI *XDS 
RELOOP 
GO  TO  40 


THE  FUNCTION  SQRT ( 1 -B/BMAX)  DOES  NOT  EXIST  FOR  B  GREATER  THAN  BMAX 
IF  PREVIOUS  STEP  IS  NOT  WITHIN  SER  OF  BMAX  INTERPOLATE  TO  FIND 
A  STEP  SIZE  THAT  WILL  GET  CLOSE  TO  BUT  NOT  EXCEED  BMAX 

70  CALL  INTERP (BB(N-2, 1) , S (N-2) , BMAX, SF, 3) 

IF (ABS (SF-S (N) ) .LT.SER)  THEN 

CALL  INTERP (BB(N-2, 1) , RSV(N-2) , BMAX, RS, 3) 

IF  (RS . LT . 1 . 03) THEN 
NMANG=IA- 1 
BMAXAN ( IA) =-1 
GOTO  218 
ENDIF 
ENDIF 

SET  UP  THE  STEP  SIZE  AND  RESET  INTEGRATION  VALUES  TO  THE  PREVIOUS 
STEP 


N=N-1 

XDS=DS 

IF (ABS (SF)  .GT. ABS  <S (N) ) )  XDS  =  0. 9* (SF-S (N) ) 

IF  THE  STEP  SIZE  IS  LESS  THAN  SER,  THE  PREVIOUS  STEP  IS  CLOSE 
ENOUGH  EXIT  TO  INVARIANT  CALCULATION 
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IF(ABS(XDS) .LT.SER)  GOTO  100 
IF (ABS (XDS) .GE.ABS(DS) ) THEN 
DS=DS/2 
ELSE 
DS=XDS 
ENDIF 

DO  80  1=1,  3 
X  < I) =XL ( I ) 

Q ( I ) =QL  ( I ) 

B { I ) =BL ( I ) 

80  CONTINUE 

85  CALL  INTGRT (X, B, BB, S,DEN) 

IF  LAST  STEP  IS  STILL  PAST  BMAX  TRY  THE  INTERPOLATION  SCHEME  AGAIN 

90  IF(BB(N, 1) .GT.BMAX)  GO  TO  70 

INTERPOLATE  TO  SEE  IF  THE  INTERPOLATION  STEP  IS  CLOSE  ENOUGH 
TO  BMAX.  IF  IT  IS  NOT,  INTERPOLATE  AGAIN  AND  TRY  TO  COME  CLOSER 

CALL  INTERP (BB(N-2, 1) , S (N-2) , BMAX, SF, 3) 

IF  WE  ARE  CLOSE  ENOUGH  EXIT  THE  INTEGRATION  LOOP 
IF (ABS (SF-S (N) ) .LT.SER)  THEN 

CALL  INTERP (BB (N-2, 1) , RSV(N-2) ,BMAX,RS, 3) 

IF  (RS.LT. 1.03) THEN 
NMANG=IA-1 
BMAXAN ( IA) =-l 
GOTO  218 
ELSE 

GOTO  100 
ENDIF 
ENDIF 

DS=DS/2 

IF (ABS (SF) . GT . ABS ( S (N) ) )  DS=CONI* (SF-S(N) ) 

CALL  INTGRT (X, B, BB, S, DEN) 

GO  TO  90 

THE  FIELD  MAXIMUM  HAS  NOW  BEEN  REACHED.  THE  STORED  VALUES 
OF  THE  MAGNETIC  FIELD  AND  THE  PATH  LENGTH  VALUES  CAN  NOW  BE 
USED  TO  EVALUATE  THE  SECOND  INVARIANT. 

100  IF (N . LT . 3 )  THEN 
XJ(IA) =0 
BMAXAN (IA) =8MAX 
GOTO  216 

ELSEIF (N.EQ. 3) THEN 
DS= . 5* (S (N-l) -S (N) ) 

CALL  INTGRT (X, B, BB, S, DEN) 

KS=2 
ELSE 
KS  =  3 
ENDIF 
NLAST=N 
DO  109  1=1,  3 
DLL (I) =BL (  I ) 

109  CONTINUE 

C  CALL  THE  ROUTINE  WHICH  DETERMINES  THE  SECOND  INVARIANT  FROM 

C  FROM  THE  STORED  VALUES 

1 1 C  CALL  INVR ( BMAX, BB, S) 

C  INTERPOLATE  GET  THE  BEST  FIT 
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CALL  INTERP (BB (N-2 , 1 ) , SXJ (N-2 ) , BMAX, XXJ, KS) 

CALL  INTERP (BB (N-2, 1 ), DEN (N-2) , BMAX, XDN, KS) 

SAVE  THE  VALUES  OF  THE  FIRST  AND  SECOND  INVARIANT 
XJ ( IA) =ABS (XXJ) 

DENSTY ( IA) =ABS (XDN) 

BMAXAN ( IA) =BMAX 


THE  INTEGRAL  HAS  NOW  BEEN  EVALUATED  FROM  THE  STARTING  POINT 
THROUGH  THE  MINIMUM  VALUE  OF  B  TO  BMAX. 

WE  MUST  INTEGRATE  THE  REST  OF  THE  LINE  -  TURN  THE  STARTING 

POINTS  AROUND  AND  RESET  THE  INITIAL  VALUES  AND  INTEGRATE  TO  THE 
OTHER  BMAX 

DEL=-DELSV 
IF ( IA . EQ . 1 ) THEN 
N=2 

BB2 (1,1) =BB (3,1) 

SXJ (1) =SXJ (3) 

S2  (1)  — S  (3) 

DEN2 ( 1 ) =DEN ( 3 ) 

DS=S (2) -S  (3) 

ELSE 

N=N2LAST 
DO  117  1=1,3 
BL (I) =BLL2  (I) 

17  CONTINUE 

DS=S2 (N) -S2 (N-l) 

CALL  INTERP (BB2 (N-2, 1) , S2 (N-2) , BMAX, SF, 3) 

IF (ABS (SF) .GT.ABS (S2 (N) ) )  THEN 
XDS=SF-S2 (N) 

IF (ABS (XDS) . LE . SER) GOTO  200 
DDS=CONI *XDS 
END  IF 
ENDIF 

CALL  STEPSZ (X2, B2, BB2 ) 

IF (ABS (BB<2, 1) -BMAX) /BMAX . LT . ERR . OR . IBEFLG.NE. 0)  GO  TO  216 
CALL  INTERP (BB (2, 1) , S (2) , BMAX, SF, 1) 

IF (ABS (SF) . LT . SER)  THEN 

CALL  INTERP (BB(2, 1) ,SXJ(2) , BMAX, XXJ, 1) 

GOTO  215 
ENDIF 

IF (ABS (SF) .LT.ABS(DS) )DS=.7*SF 
DO  120  1=1,3 
Qd)  -0. 

120  CONTINUE 
N=N2LAST 

GO  TO  140 

BEGIN  INTEGRATING  THE  SECOND  PART 
130  CALL  STEPSZ (X2, B2, BB2) 

140  CONTINUE 

CALL  INTGRT (X2,B2,BB2,S2, DEN 2 ) 

DDS=100 . 

STOP  INTEGRATION  IF  BMAX  HAS  BEEN  PASSED 
IF(BB2(N, 1) .GE.BMAX)  GO  TO  150 
IF ( IERFLG . NE . 0 ) THEN 
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IF ( IERFLG . GT . 0) THEN 
BMAXAN { I A ) =100 
ELSE 

BMAXAN (IA)  =-l 
ENDIF 

NMANG=IA-1 
GOTO  218 
ENDIF 

CALL  INTERP  (BB2 (N-2 , 1 ) , S2 (N-2 ) , BMAX, SF,  3) 

IGNORE  INTERPOLATION  IF  RESULT  IS  NOT  MONOTONIC 
DDS=100 

IF (ABS (SF) . LE .ABS (S2 (N) ) )  GOTO  130 
XDS=SF-S2 (N) 

STOP  INTEGRATION  IF  WITHIN  SER  OF  BMAX 
IF(ABS(XDS) .LT.SER)  GO  TO  200 
DDS=CONI *XDS 
GO  TO  130 

BMAX  HAS  BEEN  PASSED,  BEGIN  INTERPOLATION  SCHEME  TO  FIND  A  POINT 
CLOSE  TO  BMAX  BUT  LESS  THAN  IT. 

150  CALL  INTERP  (BB2  (N-2,  1)  ,  S2  (N-2)  ,  BMAX,  SF,  3) 

IF (ABS (SF-S2 (N) ) . LE . SER)  THEN 

CALL  INTERP (BB2 (N-2, 1) ,RSV(N-2) , BMAX, RS, 3) 

IF  (RS . LT . 1 . 03 ) THEN 
NMANG=IA-1 
BMAXAN (IA) =-l 
GOTO  218 
ENDIF 
ENDIF 

N=N- 1 
XDS=DS 

IF (ABS (SF) .GT . ABS (S2 (N) ) )  XDS=0 . 9* (SF-S2 (N) ) 

IF (ABS (XDS) .LT.SER)  GOTO  200 
IF (ABS (XDS) .GE.ABS(DS) ) THEN 
DS=DS/2 
ELSE 
DS=XDS 
ENDIF 

DO  160  1=1, 3 
X2  ( I ) =XL ( I ) 

Q  ( I) =QL ( I ) 

B2  ( I ) =BL ( I ) 

160  CONTINUE 

CALL  INTGRT (X2 , B2 , BB2 , S2 , DEN2 ) 

170  IF (BB2 (N, 1 ) . GT . BMAX)  GO  TO  150 

IF  THE  POINT  IS  LESS  THAN  BMAX  MAKE  SURE  IT  IS  CLOSE  ENOUGH.  IF 
NOT,  TRY  TO  GET  CLOSER 

CALL  INTERP ( BB2 (N-2, 1) , S2 (N-2) , BMAX, SF, 3) 

IF (ABS (SF-S2 (N) ) .LT.SER)  THEN 

CALL  INTERP (BB2 (N-2, 1) , RSV(N-2)  ,  BMAX, RS, 3) 

IF  (RS.LT. 1 . 03) THEN 
NMANG=IA-1 
BMAXAN ( I A) =-l 
GOTO  218 
ELSE 

GOTO  2  0  r- 
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ENDIF 
ENDIF 
DS=DS/2 

IF (ABS (SF)  .GT.ABS  < S2 (N) ) )  DS=CONI* (SF-S2 (N) ) 

CALL  INTGRT ( X2 , B2 , BB2 , S2 , DEN2 ) 

GO  TO  170 

INTEGRAL  IS  COMPLETE  USE  STORED  VALUES  TO  GET  INVARIANT 
200  CALL  INVR (BMAX, BB2,  S2) 

CALL  INTERP (BB2 (N-2, 1) , SXJ(N-2) , BMAX,XXJ, 3) 

CALL  INTERP (BB2 (N-2, 1) , DEN2 (N-2) , BMAX, XDN, 3) 

N2LAST=N 
DO  205  1=1, 3 
BLL2 ( I ) =BL ( I ) 

205  CONTINUE 

ADD  IN  REMAINING  CONTRIBUTION  OF  SECOND  INVARIANT 

215  XJ ( IA) =XJ ( IA) +ABS (XXJ) 

DENSTY ( IA) =DENSTY ( IA) +ABS(XDN) 

216  CONTINUE 


WE  ARE  DONE  WITH  ALL  THE  INTEGRALS  -  SET  UP  ANY  ERROR  VALUES 

18  IF (NMANG.LT. NUMANG'  THEN 

DO  219  I=NMANG+1, NUMANG 
XJ ( I ) = BMAXAN (NMANG+1 ) 

BMAXAN ( I) =BMAXAN (NMANG+1) 

DENSTY (I) = BMAXAN (NMANG+1) 

19  CONTINUE 
ENDIF 

IF  INVARIANT  EXIST  CALCULATE  L'S 
DO  220  1=1, NUMANG 

IF (BMAXAN ( I ) . GT . 0 . AND . BMAXAN ( I ) . NE . 100 ) THEN 
CALL  HILTEL (BMAXAN (I) ,XJ(I)  ,  EL ( I ) ) 

ELSEIF (BMAXAN (I) .LT.O)  THEN 

EL ( I ) =-l 

ELSE 

EL ( I ) =100 
ENDIF 

220  CONTINUE 
RETURN 
END 
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SUBROUTINE  CHECK ( PANGLE , NUMANG, IERJ 
DIMENSION  PANGLE (*) 

CHECK  TO  SEE  IF  THE  PITCH  ANGLES  ARE  BETWEEN  0  AND  90  AND  THE  THE 
PITCH  ANGLE  ARRAY  IS  MONOTONICALLY  DECREASING 
IER=0 

IF  ( PANGLE ( 1) .GT. 90. .OR. PANGLE (1) .LE.O)  IER=1 
IF  (NUMANG. GT. 1)  THEN 
DO  10  1=2, NUMANG 

IF (PANGLE (I) . GT. PANGLE (1-1) )  IER=1 
IF (PANGLE (I) . GT . 90 . . OR . PANGLE ( I ) .LE.O)  IER=1 
CONTINUE 
END  IF 
RETURN 
END 
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FUNCTION  ADENS(X) 


DETERMINE  THE  AVERAGE  ATMOSPHERIC  DENSITY  AT  A  GIVEN  ALTITUDE  IN 
GM/ (CM2  *Re) /3 

DIMENSION  X (3) 

NOMINAL  VALUE  OF  F10.7 
DATA  FI 07/114./ 

CONSTANT  THAT  CONVERTS  TO  CENTIMETERS  AND  APPLIES  THE  DIVIDE  BY 
THREE  FROM  GILL'S  METHOD 
5.7339E-3=6.371E8*2.7E-ll/3 
R=SQRT (X (1) **2+X (2) **2+X(3) **2) 

IF  (R.GT.3)  THEN 
ADENS=0 
ELSE 

A=0 . 99+ . 518*SQRT (FI  07/55) 

R=(R-1) *6371 
IF  (R. LT .110)  THEN 
ADENS=0 
ELSE 

CON= (120-R) / (A*SQRT (R-103) ) 

ADENS=5 . 7339E-3*EXP (CON) 

END  IF 
END  IF 
RETURN 
END 
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SUBROUTINE  INVR (BMAX, BB, S) 

C 

C  PURPOSE 

C  TO  CALCULATE  THE  VALUE  OF  THE  SECOND  INVARIANT 

C 

C  METHOD 

C  USE  THE  VALUES  STORED  IN  THE  S  AND  BB  ARRAYS  TO  EVALUATE  THE 

C  INTEGRAL  SQRT ( 1-BB/BMAX)  ALONG  THE  FIELD  LINE.  USE  THE 

C  SAME  INTERGRATION  METHOD  (GILLS  METHOD)  USED  IN  INTEGRATING 

C  THE  FIELD  LINE 


C 

C 

C 

C 

r 

C 

C 

C 


c 

c 

c 


INPUT  --  COMMON  BLOCK  INTPAR 

N  THE  NUMBER  OF  INTEGRATION  STEPS 

BMAX  THE  VALUE  OF  THE  MAXIMUM  MAGNETIC  FIELD  (THE  POINT 
WHERE  THE  PARTICLE  FAS  ITS  MIRROR  POINT) 

BB  A  REAL  2  DIMENSIONED  ARRAY  CONTAINING  ALL  OF  THE 

MAGNETIC  FIELD  MAGNITUDES  CALCULATED  IN  THE  FIELD  LINE 
INTEGRATION 

S  AN  ARRAY  THAT  HOLDS  THE  TOTAL  INTEGRATED  PATH  LENGHT  ALONG 

FIELD  LINE 

OUTPUT  --  COMMON  BLOCK  INTPAR 

SXJ  THE  VALUES  OF  THE  SECOND  INVARIANT  INTEGRATION  AT 
EACH  INTEGRATION  STEP.  SXJ (N)  CONTAINS  THE  BEST 
APPROXIMATION  TO  THE  VALUE  OF  THE  SECOND  INVARIANT. 

THE  SAVING  OF  THE  STEPS  PERMITS  THE  USE  OF  INTERPOLATION 
SCHEMES  TO  OBTAIN  A  MORE  ACCURATE  VALUE  OF  THE  INVARIANT 


C  CALLING  SUBROUTINES 

C  SUBROUTINE  INVARM 

C 

C  CONSTANTS 

C  P29  1 . 0-SQRT ( . 5 ) 

C  OP7  1 . 0+SQRT ( . 5) 

C 

COMMON  /INTPAR/DS,  DEL,  N,  IERFLG,  XL  ( 3)  ,  XSVUOO,  3,4), 

*RSV(100) , RMIN ( 3) , RMAG, IDSW, 

*QL ( 3 ) , Q ( 3 ) , BL ( 3 ) , SXJ (100) ,DDS 
DIMENSION  BB ( 1 00, 4 ) , S  ( 100 ) 

DIMENSION  CON (4) 

DATA  (CON (I) , 1  =  1,  4) / . 5,  .2  92  89322,  1.7  071067  8,  .5/ 

SXJ (2) =0 . 

START  THE  INTEGRATION  LOOP. 

THIS  IS  GILLS  METHOD  MADE  SIMPLE  IF  ALL  THE  POINTS  ARE  GIVEN 
CUMULATIVE  ROUND  OFF  ERROR  CONTROL  IS  NOT  IMPLEMENTED 
NN=N- 1 

DO  210  K=2 , NN 
TEMP  1  =  0 
DO  110  1=1,4 

IF (BB(K, I) .GE.BMAX)  GO  TO  110 
ROOT=SQRT (1 . -BB (K, I) /BMAX) 

TEMP  1 =TEMP 1 +CON ( I ) *ROOT 
CONTINUE 

DELS= (S (K+l) -S  (K)  ) / 3 . 

SXJ ( K+ 1 ) =SXJ (K) +TEMP 1 *DELS 
210  CONTINUE 
RETURN 
END 
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SUBROUTINE  STEPSZ (X, B, BB) 


C 

C  PURPOSE 

C  DETERMINE  THE  STEP  SIZE  FOR  THE  NEXT  INTEGRATRION  STEP 

C 

C  METHOD 

C  THE  STEP  SIZE  OF  THE  RUNGE  KUTTA  INTEGRATION  IS  A  FUNCTION 

C  OF  THE  ERROR  LIMITS,  THE  CURVATURE  OF  THE  FIELD  LINE,  THE 

C  GRADIENT  IN  THE  MAGNETIC  FIELD,  AND  THE  ESTIMATED  DISTANCE 

C  TO  THE  END  OF  THE  INTEGRATION. 

C 

C  INPUT 

C  DEL 

C 
C 

C  B 

C 

C  BL 

C 

C  BB 

C 
C 
C 

c 

C  DDS 

C 
C 

C  INPUT/OUTPUT  —  COMMON  BLOCK  INTPAR 

C  DS  ON  ENTRY  TO  THE  ROUTINE  DS  CONTAINS  THE  SIZE  OF  THE 

C  LAST  STEP.  THE  ROUTINE  RESETS  THE  VALUE  TO  THE  BEST 

C  STEP  SIZE  FOR  THE  NEXT  INTEGRATION  STEP. 

C 

C  CALLING  SUBROUTINES 

C  INVARM 

C 

C  TEMPORARY  VARIABLES 

C  CURVMN  THE  MINIMUM  ACCEPABLE  CURVATURE.  THIS  LIMITS  THE  STEP 

C  SIZE  IN  THE  VICINITY  OF  THE  EARTH  WHERE  THE  FIELD 

C  CHANGES  RAPIDLY 

C  CURV  THE  CURVATURE  OF  THE  FIELD  LINE 

C 

COMMON  / INTPAR/DS, DEL, N, IERFLG, XL < 3) , XSV ( 100 , 3,4)  , 

*RSV (100) , RMIN (3) , RMAG, IDSW, 

*QL (3) , Q (3) , BL (3) ,SXJ(100) , DDS 

DIMENSION  BB (100, 4) ,BB2  <100, 4) ,B(3),B2(3),X(3),X2(3),S(3),S2(3) 


DETERMINE  THE  MINIMUM  CURVATURE 

CURVMN=1 . 6667/  (X ( 1 ) **2+X (2) **2+X ( 3) * *2) ** ( . 75) 

DETERMINE  THE  CURVATURE  OF  THE  FIELD  BY  USING  THE  RATE  OF  CHANGE 
OF  THE  UNIT  FIELD  VECTOR  OVER  THE  LAST  STEP 

CURV=0 . 

DO  50  1=1,3 

CURV=CURV+ ( ( B ( I ) /BB(N,  1) -BL(I)  /BB(N-1,  1)  )  /DS) **2 
50  CONTINUE 

CURV=SQRT (CURV) 

CURV=AMAX1 (CURV,  CURVMN) 

SET  UP  THE  NEW  STEP  SIZE  AND  LIMIT  THE  STEP  SIZE  TO  LESS  THAN  2.8 
EARTH  RADDII  TO  PREVENT  THE  INTEGRATION  FROM  STEPPING  OUT  OF  THE 
VALID  FIELD  REGION 


—  COMMON  BLOCK  INTPAR 

A  PARAMETER  SET  UP  BY  THE  CALLING  PROGRAM  TO  SCALE  THE 
STEP  SIZE.  IT  DEPENDS  ON  THE  ERROR  LIMITS  OF  THE 
INTEGRATION. 

A  REAL  ARRAY  WHICH  CONTAINS  THE  MAGNETIC  FIELD  VECTOR 
AT  THE  CURRENT  STEP 

A  REAL  ARRAY  WHICH  CONTAINS  THE  MAGNETIC  FIELD  VECTOR 
AT  THE  PREVIOUS  STEP 
A  2  DIMENSIONED  REAL  ARRAY 

BB (N, 1)  IS  THE  MAGNETIC  FIELD  MAGNITUDE  AT  THE  CURRENT 
STEP 

BB (N-l , 1 )  IS  THE  MAGNETIC  FIELD  MAGNITUDE  AT  THE 
PREVIOUS  STEP 

THE  ESTIMATED  STEP  SIZE  REQUIRED  TO  COMPLETE  THE 
INTEGRATION 
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DS=DEL/CURV 

DS=SIGN (AMIN1 (ABS (DS) , 1 . 0) ,DS) 

IF (N . LE . 3)  DS=DS* (N*2-3) * .2 

IF  THE  DISTANCE  TO  THE  END  OF  THE  INTEGRATION  IS  SMALLER  THAN  THE 
NEW  STEP  SIZE,  SET  THE  STEP  SIZE  TO  THE  SMALLER  VALUE. 

IF (ABS (DDS) . LT .ABS (DS) )  DS=DDS 

RETURN 

END 
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c 
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c 
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c 

c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  INTGRT (X, B, BB, S, DEN) 

PURPOSE 

THIS  SUB  MODULE  PERFORMS  A  SINGLE  RUNGE-KUTTA  INTEGRATION 
STEP  AND  UPDATES  ALL  OF  THE  VARIABLES  IN  THE  INTEGRATION  LOOP 

METHOD 

PERFORM  A  SINGLE  FOURTH  ORDER  INTEGRATION  STEP  USING  GILLS 
METHOD  OF  INTEGRATION  (REF.  S.  GILL  CAMBRIDGE  PHILOSOPHICAL 
SOCIETY  PROCEEDINGS  VOL.  47,  1951) 

INPUT  —  COMMON  BLOCK  INTPAR 

DS  THE  INTEGRATION  STEP  SIZE  IN  UNITS  OF  EARTH  RADII. 

THE  INTEGRATION  MOVES  THE  SPACE  COORDINATE  A  DISTANCE 
DS  ALONG  THE  MAGNETIC  FIELD  LINE.  IF  DS  IS  POSITIVE, 
MOTION  IS  IN  THE  DIRECTION  OF  THE  FIELD.  IF  DS  IS 
NEGATIVE  MOTION  IS  ANTI-PARALLEL  TO  THE  FIELD. 

INPUT/OUTPUT  --  COMMON  BLOCK  INTPAR 

N  THE  INTEGRATION  STEP  NUMBER.  IT  IS  INCREMENTED  BY 

ONE  AT  THE  END  OF  THIS  ROUTINE.  (NOTE  N=2  IS  THE 
BEGINNING  OF  THE  INTEGRATION) 

X  A  REAL  ARRAY  GIVING  THE  VECTOR  LOCATION  OF  THE 

INTEGRATION  VARIABLE. 

INPUT  -  THE  INITIAL  POSITION  PRIOR  TO  THE  INTEGRATION 
STEP 

OUTPUT-  THE  FINAL  VALUE  AFTER  THE  INTEGRATION  STEP 
B  A  REAL  ARRAY  CONTAINING  THE  VECTOR  MAGNETIC  FIELD 

IN  GAUSS 

INPUT  -  THE  VECTOR  FIELD  BEFOR  THE  INTEGRATION  STEP 
OUTPUT-  THE  VECTOR  FIELD  AFTER  THE  STEP 
Q  A  REAL  ARRAY  CONTAINING  AN  ERROR  CONTROL  VARIABLE 

USED  BY  GILLS  INTEGRATION  METHOD 
INPUT  -  ERROR  FROM  PREVIOUS  STEP 

OUTPUT-  ERROR  AFTER  PRESENT  STEP  FOR  INPUT  TO  SUBSEQUENT 
STEPS 

OUTPUT  --  COMMON  BLOCK  INPTAR 

S  A  REAL  ARRAY  WHICH  SAVES  EACH  OF  THE  DISTANCES  (SINCE 

THE  START  OF  THE  INTEGRATION)  ALONG  THE  MAGNETIC  FIELD 
LINE. 

S (2) =0 

S (N+l) =S (N) +DS  ETC. 

XSV  A  REAL  3  DIMENSIONED  ARRAY  WHICH  SAVES  THE  VECTOR 

POSITION  IN  EARTH  RADII  FOR  EACH  OF  THE  INTEGRATION 
STEPS.  XSV(N, 1,1),  XSV (N, 2 , 1 ) ,  XSV(N,3,1)  ARE  VECTOR 
CARTESIAN  POSITION  COORDINATES  CORESPONDING  TO  POSITION 
S (N)  ON  THE  FIELD  LINE 

BB  A  REAL  2  DIMENSIONED  ARRAY  WHCIH  SAVES  THE  MAGNITUDE 

OF  THE  MAGNETIC  FIELD  FORM  EACH  INTEGRATION  STEP. 

BB (N, 1 )  IS  MAGNETIC  FIELD  VALUE  AT  DISTANCE  S (N) . 

BB (N-l , 2 ) ,  BB (N-l, 3) ,  BB(N-1,4)  ARE  THE  INTERMEDIATE 
VALUES  OF  THE  FIELD  USED  BY  GILLS  METHOD  TO  GET  FROM 
BB (N-l,  1)  TO  BB (N, 1)  . 

XL  A  REAL  ARRAY  WHICH  SAVES  THE  INITIAL  POSITION  VALUES 

PRIOR  TO  STARTING  THE  INTEGRATION  STEP 
BL  A  REAL  ARRAY  WHICH  SAVES  THE  VECTOR  MAGNETIC  FIELD 

VALUES  PRIOR  TO  STARTING  THE  INTEGRATION  STEP 
QL  A  REAL  ARRAY  WHICH  SAVES  THE  INITIAL  VALUES  OF  THE 

ERROR  CONTROL  VARIABLE 

IERFLG  AN  ERROR  CONTROL  INDICATOR  WHICH  IS  USED  BY  THE  CALLING 
PROGRAM  TO  CONTROL  THE  PROGRAM  FLOW 
IERFLG  =  0  NO  ERROR 

IERFLG  =  1  INTEGRATION  IS  OUTSIDE  VALID  FIELD  LIMITS 
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OR  THE  MAXIMUM  STEP  NUMBER  (100)  HAS  BEEN 
REACHED . 

CONSTANTS 

P29  1  .  0-SQRT (0.5) 

OP  7  1  .OfSQRT(O.S) 

VARIABLES 

P5DS  .5  *  STEP  SIZE 
P29DS  ( 1. O-SQRT  (0. 5) )  *  STEP  SIZE 
OP7DS  ( 1 . 0-SQRT (0 . 5) )  *  STEP  SIZE 

RR, SS  REAL  ARRAYS  USED  BY  GILLS  METHOD  TO  MINIMIZE  COMPUTER 
TIME  AND  MINIMIZE  ROUNDOFF  ERROR 

CALLING  SUBROUTINES 
SUBROUTINE  INVARM 

SUBROUTINES  REQUIRED 
SUBROUTINE  BMNEXT 

COMMON/BXYZCM/ YEAR, DAYYR, UT, KODE,  JSW 

COMMON  / INTPAR/DS, DEL, N, IERFLG, XL ( 3)  ,  XSV(100,  3,4), 

*RSV(100) , RM I N ( 3 ) , RMAG, IDSW, 

*QL (3) , Q (3) , BL ( 3) , SXJ(IOO) ,DDS 
DIMENSION  BB(100, 4) ,B(3) , X ( 3 ) ,S(100) , DEN  (100) 

DIMENSION  SS ( 3 ) ,RR(3) 

DATA  P2 9, OP  7/ .29289322,  1 .70710678/ 

IERFLG=0 

SAVE  THE  INITIAL  VALUES.  THESE  INITIAL  VALUES  MAY  BE  NEEDED  IF 
IF  THE  INTEGRATION  STEP  IS  UNSUCCESSFUL  (GOES  TOO  FAR)  AND  THE 
STEP  MUST  BE  REPEATED. 

DO  65  1  =  1,  3 
XL ( I ) =X ( I ) 

QL ( I ) =Q ( I ) 

BL ( I ) =B ( I ) 

63  CONTINUE 

SET  UP  THE  CONSTANST  NEEDED  BY  THE  INTEGRATION  LOOP 

P  5DS= . 5  *DS 
P29DS=P29*D3 
OP7DS  =  OP7  *CS 

BEGIN  GILLS  METHOD  (GILL  1951)  OF  FOURTH  ORDER  INTEGRATION 

TEMP2=P5DS  *  ADEN'S  (X) 

DO  70  1=1,3 

SS  ( I)  =P  5DS  *  B  (I)  /BB  (N,  1) 

RR (I) =SS (I) -Q(I) 

X ( I ) =X ( I ) +RR  ( I ) 

Q(I) =Q(I) +3. ‘RR (I) -SS(I) 

XSV (N, I , 2 ) =X ( I ) 

"0  CONTINUE 

TEMP2  =  TEMP2  +  P2 9DS ‘ADEN’S  ( X) 

CALL  BMNEXT (X, B, BB(N,  2)  ) 

DO  71  1=1,3 

SS ( I ) =P2  9DS  *B ( I ) /BB (N, 2 ) 

RR ( I) =SS (I) -?2  9*Q ( I ) 

x ( I ) =x  ( i )  ‘F p  ( : ) 

Q  (I)  =Q  (I)  +3  .  *  RR  ( I )  -SS  (1) 

XSV (N, I, 3) — X(I) 

71  CONTINUE 

TEMP2  =  TEMP2  1  "  P  70S  *  ADEN  S  (X) 
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CALL  BMNEXT (X,  B,  BB  (N,  3)  ) 

DO  72  1=1,3 

SS;i)=0P7DS*B(I) /BB(N,3) 

RR  < I ) =SS  ( I ) -OP 7  *Q  ( I ) 

X(I)=X(I) +RR ( I ) 

Q (I) =Q ( I) +3 . *RR (I) -SS (I) 

XSV  (N,  I,  4)  =X  (I) 

72  CONTINUE 

TEMP2=TEMP2  +P5DS  * ADENS (X) 

CALL  BMNEXT (X,B,BB(N,  4)  ) 

DO  73  1=1,3 

SS ( I) =P5DS*B ( I ) /BB (N, 4) 

RR  ( I )  =  ( SS ( I ) -Q ( I ) ) / 3  . 

X ( I ) =X ( I ) +RR ( I ) 

Q(I)=Q(I)+3.*RR(I)-SS(I) 

XSV  (N+l,  I,  1)  =X  (I) 

73  CONTINUE 
N=N+1 

SAVE  THE  CURRENT  DISTANCE  ALONG  THE  FIELD  LINE 

S (N) =S ( N— 1 ) +DS 
DEN (N) =DEN (N-l ) +TEMP2 

OBTAIN  THE  CURRENT  VALUES  OF  THE  MAGNETIC  FIELD 

CALL  BMNEXT (X,B,BB(N,  1) ) 

IF  N  IS  TOO  BIG,  SET  ERROR  FLAG 
IF(N.GE.IOO)  IERFLG=1 

IF  OUTSIDE  INTEGRATION  LIMITS  SET  ERROR  FLAG 
R=X ( 1 ) **2+X (2) **2+X (3) **2 
RSV (N) =SQRT (R) 

C  IF  EXTERNAL  FIELD  IS  USED  STAY  WITHIN  VALID  REGION 

IF( <R.GT. 144. .OR. BB(N, 1) . LT . 0.00015) . AND . JSW . GE . 0 ) 
C  IF  BELOW  EARTHS  SURFACE  SET  FLAG  NEGATIVE 

IF (RSV (N) .LT.1.03)  IERFLG=-1 
RETURN 
END 


IERFLG=1 
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SUBROUTINE  INTERP (BB, CC, D, E, J) 

C 

C  PURPOSE 

C  INTERPOLATION  ROUTINE 

C 

C  METHOD 

C  GIVEN  A  SET  OF  THREE  X,  Y  POINT  PAIRS,  INTERP  FINDS  THE  SOLUTION 

C  TO  THE  THREE  LINEAR  EQUATIONS  EXPRESSING  THE  LOGARITHM  OF  THE 

C  DEPENDENT  VARIABLE  Y  AS  A  SECOND  ORDER  POLYNOMIAL  OF  THE 

C  INDEPENDENT  VARIABLE  X.  (LOG  Y  =  A*X**2  +  B*X  +C) 

C  USING  THE  BINOMIAL  FORMULA,  X  CAN  THEN  BE  EVALUATED  AT  A 

C  SPECIFIED  VALUE  OF  Y1 

C  X  =  (-B  +  -  SQRT (B**2-4*A* (C-LOG(Yl)) ) ) / (2*A) 

ARGUMENT  LIST 

A  REAL  ARRAY  CONTAINING  THE  THREE  VALUES  OF  THE 
DEPENDENT  VARIABLE 

A  REAL  ARRAY  CONTAINING  THE  THREE  CORESPONDING  VALUES 
OF  THE  INDEPENDENT  VARIABLE 
A  FLOW  CONTROL  VARIABLE 
IF  J  IS  LESS  THAN  0 

FIT  THE  POLYNOMIAL  TO  CC  AND  BB  AND  FIND  THE  MINIMUM 
VALUE  OF  THE  DEPENDENT  VARIABLE 
IF  J  IS  GRATER  THAN  0 

USE  THE  BINOMIAL  FORMULA  TO  TO  FIND  THE  VALUE  OF 
THE  INDEPENDENT  VARIABLE  WHEN  THE  DEPENDENT  VARIABLE 
HAS  THE  VALUE  D. CHOOSE  THE  ROOT  THAT  IS  CLOSEST  TO  CC(J) 
WHEN  J  IS  GREATER  THAN  ZERO,  D  IS  USED  FOR  INPUT. 

IT  IS  THE  VALUE  OF  THE  DEPENDENT  VARIABLE  WHERE  THE 
SOLUTION  TO  THE  DEPENDENT  VARIABLE  IS  WANTED 


C 

C 

c 


INPUT 

BB 

CC 

J 


C  OUTPUT  --  ARGUMENT  LIST 

C  D  WHEN  J  IS  LESS  THAN  0,  D  OUTPUTS  THE  VALUE  OF  THE 

C  DEPENDENT  VARIABLE  WHERE  THE  FUNCTION  IS  A  MINIMUM 

C  E  WHEN  J  IS  LESS  THAN  0,  E  OUTPUTS  THE  VALUE  OF  THE 

C  INDEPENDENT  VARIABLE  WHERE  THE  FUNCTION  IS  A  MINIMUM 

C  WHEN  J  IS  GREATER  THAN  0,  E  OUPUTS  THE  VALUE  OF  THE 

C  INDEPENDENT  VARIABLE  WHERE  THE  FUNCTION  HAS  THE  VALUE  D 


CALLING  SUBROUTINES 
SUBROUTINE  INVARM 


VARIABLES 

X2 , X 3 , Y 1 , Y 2 , Y 3 , DD  ARE  USED 
TO  MINIMIZE  COMPUTER 
A, B, C  THE  THREE  POLYNOMIAL 
DIS  B*'2-4‘A*C 
SA, SB  THE  TWO  ROOTS  OF  THE 


BY  THE  LINEAR  EQUATION  SOLUTION 
TIME 

COOEFICIENTS 

POLYNOMIAL 


DIMENSION  BB ( 3 )  ,  CC ( 3 ) 

REAL* 6  Y1,Y2,Y3,X2,X3, DD, A, B, C,  DIS 


SET  UP  THE  INITIAL  VARIABLES,  MOVE  THE  ORIGIN  OF  THE  INDEPENDENT 
VARIABLE  TO  CC ( I ) 


Yl=ALOG (BB ( 1 ) ) 

Y2=ALOG (BB (2) ) 

Y3=ALOG (BB ( 3) ) 

X2=CC (2) -CC  (I ) 

X3=CC ( 3 ) -CC(1) 

SOLVE  THE  LINEAR  EQUATIONS 
DD=  (X3-X2 )  *X2  *X3 
IF  (DD . EQ . 0 ) 7"’  N 
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nooono  on  on  no  non  ooo 


IF(J.LT.O)  THEN 
D=BB (2 ) 

ELSE 

E=CC ( J) 

ENDIF 

RETURN 

ENDIF 

A= (X3* ( Y1-Y2 ) +X2  * (Y3-Y1) ) /DD 
B« (X3**2* (Y2-Y1) -X2**2* (Y3-Y1) ) /DD 

IF  J  THE  FLOW  CONTROL  VARIABLE  IS  LESS  THAN  ZERO  BRANCH  TO 
MINIMUM  EVALUATION  ROUTINE 
IF ( J . LT . 0) GO  TO  100 
C=Yl-DLOG (D) 

DIS=rB**2-4 .  *A*C 

IF  DIS  IS  NEGATIVE  NO  SOLUTION  EXIST,  EXCHANGE  DEPENDENT  AND 
INDEPENDENT  VARIABLE  ROLES  AND  TRY  ANOTHER  SOLUTION 
IF (DIS . LE . 0 . )  GO  TO  2C0 
DIS=DSQRT (DIS) 

OBTAIN  THE  TWO  ROOTS 
SA= (-B+DIS) / (2 . *A) +CC(1) 

SB= (-B-DIS) / (2 . * A) +CC (1) 

E=SA 


FIND  THE  ROOT  CLOSEST  TO  CC(J) 

IF (ABS (SB-CC (J) ) . LT . ABS (SA-CC (J) ) )  E=SB 
RETURN 

FIND  THE  VALUES  AT  THE  MINIMUM 
100  X— B/(2.*A) 

E=X+CC  (1) 

XM=A*X*  *2+B*X+Yl 
D=EXP (XM) 

RETURN 

ALTERNATE  INTERPOLATION  SCHEME  PLACED  HERE  AS  A  SAFEGUARD 
AGAINST  A  STRANGE  FIELD  CONFIGURATION  CAUSING  AN  IMAGINARY 
SOLUTION  (EXCHANGE  THE  ROLES  OF  DEPENDENT  AND  INDEPENDENT 
VARIABLES) 

200  Y1=CC ( 1 ) 

Y2=CC (2) 

Y3=CC ( 3) 

X2=BB (2 ) -BB ( 1 ) 

X3=BB (3) -BB(1) 

DD= (X3-X2) *X2  *X3 
IF (DD . EQ . 0) THEN 
E=CC (J) 

RETURN 

ENDIF 

A=  (X3*  (Y1-Y2)  +-X2*  (Y3-Y1)  )  /DD 

B= ( X3  *  *2  * (Y2-Y1) -X2**2* (Y3-Y1) ) /DD 

DX=D-BB  ( 1 ) 

E= ( A*DX+B) *Dx+Yl 

RETURN 

END 
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SUBROUTINE  MG LONG (X, S, SF, XMLONG, RMIN, RMAG) 

PURPOSE 

TO  DETERMINE  THE  MAGNETIC  LONGITUDE  OF  THE  MINIMUM  B  LOCATION 
OF  THE  MAGNETIC  FIELD  LINE 

METHOD 

GIVEN  A  LOCUS  OF  POSITIONS  ALONG  A  FIELD  LINE  AS  A  FUNCTION 
OF  THE  SCALAR  DISTANCE  ALONG  THE  FIELD  LINE  AND  GIVEN  THE 
SCALAR  DISTANCE  WHERE  THE  FIELD  IS  A  MINIMUM,  THE  ROUTINE 
FINDS  THE  VECTOR  POSITION  OF  THE  MINIMUM.  IT  THEN  TRANSFORMS 
THIS  MINIMUM  TO  OFFSET  DIPOLE  COORDINATES  AND  CALCULATES 
THE  MAGNETIC  LONGITUDE  OF  THE  MINIMUM 

NOTE*  ******THE  CONSTANT  ISWTCH  IS  SET  BY  A  DATA  STATEMENT, 

IF  IT  IS  SET  TO  ZERO  XMLONG  IS  CALCULATED  USING  A  CENTERED 
DIPOLE  COORDINATE  SYSTEM  WITH  ZERO  LONGITUDE  AT  69  DEGREES 
WEST  GEOGRAPHIC.  IF  ISWTCH  IS  SET  NON-ZERO,  AN  OFFSET  DIPOLE 
COORDINATE  SYSTEM  IS  USED  WITH  XMLONG=0  GOING  THROUGH 
GREENWHICH 

INPUT  —  ARGUMENT  LIST 

X  A  REAL  2  DIMENSIONED  ARRAY  CONTAINING  THE  LOCUS  OF 

POINTS  ALONG  A  FIELD  LINE 

X(l,l),  X ( 1 , 2 )  AND  X ( 1 , 3)  ARE  THE  X,  Y,  Z  VALUES 
(RIGHT  HANDED  CARTESIAN  COORDINATES)  AT  THE  FIRST 
POINT,  X (2, 1) ,  X (2, 2)  AND  X(2,3)  THE  SECOND  LOCATION 
AND  X ( 3 , 1 ) ,  X (3, 2)  AND  X<3,3)  ARE  AT  THE  THIRD  LOCATION 
THE  FIRST  DIMENSION  OF  X  MUST  BE  THE  SAME  AS  THE 
CALLING  PROGRAMS  DIMENSION  -  IN  THIS  CASE  IT  IS  100 
S  A  REAL  ARRAY  CONTAINING  THE  SCALAR  DISTANCE  ALONG  THE 

FIELD  LINE  IN  EARTH  RADII.  S(l)  IS  THE  SCALAR  DISTANCE 
TO  THE  X(l,l),  X<1, 2),  X(l,3)  POINT  FROM  THE  START 
OF  THE  INTEGRATION,  S(2)  IS  THE  DISTANCE  TO  X(2,l),... 
SF  THE  SCALAR  DISTANCE  TO  THE  MAGNETIC  MINIMUM 

OUTPUT  --  ARGUMENT  LIST 

XMLONG  THE  MAGNETIC  LONGITUDE  (IN  DEGREES)  OF  THE  MINIMUM 
OF  THE  MAGNETIC  LINE  OF  FORCE 

IF  ISWTCH  IS  ZERO,  THE  ZERO  OF  MAGNETIC  LONGITUDE  IS 
ALONG  69  DEG  WEST  GEOGRAPHIC 

IF  ISWTCH  IS  NOT  ZERO,  THE  ZERO  OF  MAGNETIC  LONGITUDE 
IS  THROUGH  GREENWHICH 


r 

C 

V- 

C 

C 

c 


c 

c 


c 

c 

c 


CONSTANTS 

DX  THE  3  VECTOR  COMPONENTS  OF  THE  LOCATION  OF  THE  OFFSET 

DIPOLE  IN  EARTH  RADII  (GEOGRAPHIC  CARTESIAN  COORDS) 
A22-A34  TRANSFORMATION  MATRIX  TO  OFFSET  DIPOLE  COORDS.  FIRST 
ROTATE  ABOUT  THE  GEOGRAPHIC  Z  AXIS,  TO  THE  MERIDIAN 
CONTAINING  THE  OFFSET  DIPOLE,  THEN  ABOUT  THE  NEW  Y  AXIS 
TO  THE  LATITUDE  CONTAINING  THE  OFFSET  DIPOLE  AND  THEN 
ABOUT  THE  NEW  Z  AXIS  SUCH  THAT  THE  ZERO  OF  LONGITUDE 
PASSES  THROUGH  GREENWHICH 
ISWTCH  A  FLOW  CONTROL  CONSTANT 

IF  SET  TO  ZERO  BY  THE  DATA  STATEMENT  USE  CENTERED  DIPOLE 
COORDINATES 

IF  SET  NON-ZERO  USE  OFFSET  DIPOLE  COORDINATES 
SIN  D  SINE  OF  THE  COLATITUDE  OF  THE  CENTERED  DIPOLE  AXIS 
COS  D  COSINE  OF  THE  COLATITUDE  OF  THE  CENTERED  DIPOLE  AXIS 
S69  SINE  OF  69  DEGREES 

C69  COSINE  OF  69  DEGREES 

TEMPORARY  VARIABLES 

XF, XI , X2 , Y1 , Y2 , Y3 , A, B, DD  THESE  VARIABLE  ARE  USED  IN  THE 

INTERPOLATION  LOOP  TO  MINIMIZE  THE  NUMBER  OF  MEMORY 
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C  REFERENCES  AND  TO  MINIMIZE  THE  NUMBER  OF  MULTIPIES 

C  XT  A  REAL  ARRAY  HOLDING  THE  LOCATION  OF  THE  MINIMUM  AND 

C  LATER  THE  OFFSET  MINIMUM  OF  THE  FIELD  LINE 

C  XP, YP  THE  POSITION  OF  THE  MINIMUM  IN  OFFSET  MAGNETIC  CORDS. 

DIMENSION  X (100, 3) , S (100) ,DX{3) , XT (3) , RMIN (3) 

DATA  DX ( 1 ) , DX (2 ) ,  DX  ( 3) /0 . 057 6, -0 . 032 1 , -0 . 01 84 / 

DATA  A22, A23, A24, A32, A33, A34/0 . 97056, 0.23948,-0.02556, 
*-0.22969,0.95232,0.20082/ 

DATA  SIND, COSD,S69,C69/. 2027872954, .9792228106, .9335804265, 
*.3583679495/ 

C 

C 

C  ******  *  *SET  UP  THE  FLOW  CONTROL  SWITCH* *  *  *  *  *  *  *  *  * 

C  COORDINATE  SYSTEM  DEFINITION  USED.  (SEE  METHOD) 

DATA  ISWTCH/1/ 

C 

C  BEGIN  QUADRATIC  INTERPOLATION 

XF=SF-S (1) 

X2=S (2) -S ( 1 ) 

X3=S  (3) -S ( 1) 

DD=  (X3-X2) *X2*X3 
C 

C  INTERPOLATE  EACH  COMPONENT  SEPERATELY 

DO  10  1-1,3 
Y1=X  (1,1) 

Y2=X (2, I) 

Y3=X (3, I) 

A=  <X3* ( Y1-Y2 ) +X2* (Y3-Y1) ) /DD 
B= (X3**2* (Y2-Y1) -X2**2* (Y3-Y1) ) /DD 

C 

C  EVALUATE  THE  POSITION  OF  THE  MINIMUM 

XT ( I )  =  ( A*XF +B) *XF +Y1 
RMIN ( I ) =XT ( I ) 

10  CONTINUE 

RMAG2=XT ( 1 ) *  *2+XT (2 ) **2+XT<3) **2 
RMAG=SQRT (RMAG2 ) 

C 

C  IF  ISWTCH  IS  ZERO  GO  TO  CENTERED  DIPOLE  DEFINITION 

IF  ( ISWTCH. EQ. 0)  GO  TO  30 
C 

C  ADD  IN  THE  DIPOLE  OFFSET 

DO  20  1=1,3 
20  XT ( I ) =XT  < I ) +DX ( I ) 

C 

C  TRANSFORM  TO  OFFSET  DIPOLE  COORDINATES  AND  EVALUATE  THE  LONGITUDE 

XP=A22*XT (1) +A23*XT (2) +A24*XT (3) 

YP=A32*XT (1) +A33*XT (2) +A34*XT (3) 

GO  TO  40 
C 

C  TRANSFORM  TO  CENTERED  DIPOLE  COORDINATES 

30  XP= (XT (1) *C69-XT (2) *S69) *COSD-XT (3) *SIND 
YP=XT ( 1 ) *S69+XT (2) *C69 
C 

C  CALCULATE  MAGNETIC  LONGITUDE 

40  XMLONG=ATAN2 (YP, XP) *57 .2957795 

IF (XMLONG . LT . 0 . )  XMLONG=XMLONG+3 60 . 

RETURN 

END 
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SUBROUTINE  HILTEL  (B,XI,VL) 

C  PURPOSE 

C  CALCULATE  THE  L  VALUE 

C  THE  ORIGINAL  MCILWAIN  L  EXPANSION  GIVEN  BY  THE  OLD 

C  SUBROUTINE  CARMEL  HAS  BEEN  REPLACED  BY  HILTONS  SIMPLER 

C  EXPANSION.  DIFFERENCES  BETWEEN  HILTONS  AND  MCILWAINS 

C  EXPANSION  ARE  TYPICALLY  LESS  THAN  .01  PERCENT. 

C 

C  METHOD 

C  SEE  J.  HILTON,  J.  GEOPHYS.  RES.  76,  6952  (1971) 

C 

C  INPUT  --  CALLING  SEQUENCE 

C  B  THE  MAGNETIC  FIELD  AT  THE  PARTICLE  MIRROR  POINT 

C  XI  THE  SECOND  INVARIANT  EVALUATED  BETWEEN  MIRROR  POINTS 

C  EXPRESSED  IN  UNITS  OF  EARTH  RADII 

C 

C  OUTPUT  --  CALLING  SEQUENCE 

C  VL  THE  L  VALUE 

C 

C  THE  NEXT  STATEMENT  CONTAINS  THE  ORIGINAL  McILWAIN  MOMENT 

C  DATA  XM/. 311653/ 

C  USE  THE  DIPOLE  MOMENT  CALCULATED  FROM  THE  CURRENT  FIELD  MODEL 

COMMON  /MOMENT /XM 
IF (XI .GT . 1 . OE-36)  GO  TO  10 
VL= (XM/B) **  (1./3.) 

RETURN 

1C  X=XI* (B/XM) **  (1./3.) 

V=1 ,+X* (1 . 35047+X* ( .465376+. 047 54 55 *X) ) 

VL= (V*XM/B)  **  ( 1 . / 3  .  ) 

C  END  COMPUTE  L 

RETURN 
END 
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